About the code

Experiment: FSAReward (Ivan Grahek, Antonio Schettino, Gilles Pourtois, Ernst Koster, & Søren Andersen) (*: co-first authors) Code written by: Ivan Grahek & Antonio Schettino (2016-2018) Description: Code for the analysis of EEG data for Experiment 1 of the SSVEP - reward project.

This is an R Markdown document. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.

EEG

Topography & spectra

Topography and spectra

Import data

# # Clear environemnt and import data------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

# clear the environment
rm(list=ls())
# clear the console
cat("\014")

#load packages and install them if they're not installed
if (!require("pacman")) install.packages("pacman")
pacman::p_load(reshape2,yarrr,BayesFactor,plyr,ez,schoRsch,brms,knitr,broom, brmstools, BEST, tidyverse, here, HDInterval)
# set seed
set.seed(42)
# Set working directory
#setwd(here())
# import data
data.raw = read.csv(file = here("EEG_preprocessing/movement","grandAverage_amplitudes.csv"),header=TRUE,na.strings="NaN") 

# Rename columns
colnames(data.raw) = c("Subject","Frequency","BslnRedAttendedNomov","BslnBlueAttendedNomov","AcqRedAttendedNomov","AcqBlueAttendedNomov","ExtRedAttendedNomov","ExtBlueAttendedNomov","BslnRedAttendedMov","BslnBlueAttendedMov","AcqRedAttendedMov","AcqBlueAttendedMov","ExtRedAttendedMov","ExtBlueAttendedMov")

# # Reshape to long format - take only no movement trials
data = melt(data.raw,id.vars=c("Subject","Frequency"),
            measure.vars=c("BslnRedAttendedNomov","BslnBlueAttendedNomov","AcqRedAttendedNomov","AcqBlueAttendedNomov","ExtRedAttendedNomov","ExtBlueAttendedNomov","BslnRedAttendedMov","BslnBlueAttendedMov","AcqRedAttendedMov","AcqBlueAttendedMov","ExtRedAttendedMov","ExtBlueAttendedMov"),
            variable.name="Condition",value.name="Amplitude")

# Sort the new dataframe by participant name
data = data[order(data$Subject),]

# Split the variable Condition based on capital letters
data$Condition = gsub("(?!^)(?=[[:upper:]])", " ", data$Condition, perl=T)

# Split the variable condition into multiple variables
Conditions = colsplit(data$Condition, pattern="\\s+",names = c('ExpPhase', 'ColorMoved',"attended","no","moved"))

# Add the variable defining which color is rewarded based on the participant number
data$RewardedColor = ifelse(data$Subject%%2==0,"Blue","Red") # if participant number is even, blue was rewarded

# Add the Conditions needed to the dataset
data$ExpPhase = Conditions[,1]
data$AttendedColor = Conditions[,2]
data$Movement = Conditions[,4]

# Switch the Frequency to the color
data$RecordedFrequency = ifelse(data$Frequency==10,"Blue","Red") # if the recorded frequency is 10Hz assign Blue (color flickering at 10Hz), otherwise assign Red (color flickering at 12Hz)

# Make a new condition based on the attended color and the rewarded color
data$Condition = ifelse(data$AttendedColor==data$RewardedColor, "High_Rew","Low_Rew")

# Make a new condition based on the attended color and the recorded frequency
data$Attention = ifelse(data$AttendedColor==data$RecordedFrequency, "Att","NotAtt")

# Make a new condition based the Condition and the Attention
data$RecordingAndCondition = with(data, paste0(Condition,"_",Attention))

# Select variables which we want to keep
data = subset(data, select=c("Subject","RewardedColor","ExpPhase","AttendedColor","Condition","RecordedFrequency","Attention","RecordingAndCondition","Movement","Amplitude"))

# Sort the data 
data = data[with(data, order(Subject)), ]

# Make a new variable with mean amplitude across all conditions for each participant and each frequency (normalization)
data = ddply(data,.(Subject,RecordedFrequency),transform,
             MeanAmplitude = mean(Amplitude[ExpPhase=="Bsln"],na.rm=TRUE),
             SDAmplitude =   sd(Amplitude,na.rm=TRUE))

# Divide amplitudes in each Subject, Frequency, and Condition by the Mean Amplitude
data$Amplitude = data$Amplitude/data$MeanAmplitude

# Calculate the attention indexes - Selectivity (attended-unattended) & total enhancement (attended+unattended) (Andersen & Muller, 2010, PNAS)
data.diff = ddply(data, .(Subject,ExpPhase,Condition), transform, Selectivity = Amplitude[Attention=="Att"]-Amplitude[Attention=="NotAtt"],TotalEnhancement=Amplitude[Attention=="Att"]+Amplitude[Attention=="NotAtt"])
# Delete the Attention column and rows which are not necessary (indexes repeated twice)
data.diff = subset(data.diff,Attention=="Att") #keep only Att as it is equal to NotAtt
data.diff$Attention = NULL  #drop the Attention column

# Sort the data 
data.diff$ExpPhase = factor(data.diff$ExpPhase, levels = c("Bsln","Acq","Ext"))
data.diff = data.diff[order(data.diff$Subject,data.diff$Condition,data.diff$ExpPhase),]

# Calculate the reward index - High reward minus Low reward
data.reward = ddply(data, .(Subject,ExpPhase,Attention), transform, Reward = Amplitude[Condition=="High_Rew"]-Amplitude[Condition=="Low_Rew"])
# Delete the Attention column and rows which are not necessary (indexes repeated twice)
data.reward = subset(data.reward,Condition=="High_Rew") #keep only Att as it is equal to NotAtt
data.reward$Condition = NULL  #drop the Condition column

# Sort the data 
data.reward$ExpPhase = factor(data.reward$ExpPhase, levels = c("Bsln","Acq","Ext"))
data.reward = data.reward[order(data.reward$Subject,data.reward$Attention,data.reward$ExpPhase),]

# Convert variables to be used in analyses into factors
data[c("Subject", "Condition","ExpPhase", "RewardedColor", "Attention", "RecordingAndCondition")] = 
  lapply(data[c("Subject", "Condition","ExpPhase", "RewardedColor", "Attention", "RecordingAndCondition")], factor)

data.diff[c("Subject", "Condition","ExpPhase", "RecordingAndCondition")] = 
  lapply(data.diff[c("Subject", "Condition","ExpPhase",  "RecordingAndCondition")], factor)

# Save the final data into a new variable
data_all = data
data_all = subset(data_all, Subject!=14)

All trials

Means - raw data

summary = ddply(data,.(Attention,ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(Amplitude, na.rm = TRUE), digits = 2), " [", round(hdi(Amplitude)[[1]], digits = 2), " ", round(hdi(Amplitude)[[2]], digits = 2), "]")))

names(summary) = c("Attention", "Reward phase", "Reward probability", "Amplitude")

summary$Attention = recode(summary$Attention,
                           "Att" = "Attended",
                           "NotAtt" = "Unattended")

summary$`Reward phase` = recode(summary$`Reward phase`,
                           "Acq" = "Acquisition",
                           "Bsln" = "Baseline",
                           "Ext" = "Extinction")

summary$`Reward probability` = recode(summary$`Reward probability`,
                           "High_Rew" = "High",
                           "Low_Rew" = "Low")

summary = as.data.frame(summary)

summary$`Reward phase` = factor(summary$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
summary = summary[order(summary$Attention,summary$`Reward phase`,summary$`Reward probability`),]
row.names(summary) = NULL

kable(summary, caption = "Amplitudes per condition")
Amplitudes per condition
Attention Reward phase Reward probability Amplitude
Attended Baseline High 1.13 [ 0.92 1.52 ]
Attended Baseline Low 1.13 [ 0.86 1.52 ]
Attended Acquisition High 1.16 [ 0.8 1.6 ]
Attended Acquisition Low 1.13 [ 0.76 1.71 ]
Attended Extinction High 1.13 [ 0.61 1.61 ]
Attended Extinction Low 1.13 [ 0.59 1.84 ]
Unattended Baseline High 0.87 [ 0.47 1.17 ]
Unattended Baseline Low 0.87 [ 0.49 1.11 ]
Unattended Acquisition High 0.91 [ 0.54 1.38 ]
Unattended Acquisition Low 0.89 [ 0.5 1.28 ]
Unattended Extinction High 0.88 [ 0.48 1.23 ]
Unattended Extinction Low 0.91 [ 0.44 1.42 ]

Plots - raw data

# Plot amplitude across experiment phases------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# prepare data for plotting
dataPlot = data

# rename variables
colnames(dataPlot)[colnames(dataPlot)=="ExpPhase"] <- "Reward phase"
colnames(dataPlot)[colnames(dataPlot)=="Condition"] <- "Reward probability"

# rename conditions
dataPlot$`Reward phase` = recode(dataPlot$`Reward phase`,
                                  "Acq" = "Acquisition",
                                  "Bsln" = "Baseline",
                                  "Ext" = "Extinction")

dataPlot$`Reward probability` = recode(dataPlot$`Reward probability`,
                                        "High_Rew" = "High",
                                        "Low_Rew" = "Low")

#order
dataPlot$`Reward phase` = factor(dataPlot$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward phase`,dataPlot$`Reward probability`),]


plottingConditions = c("Attended","Unattended" )
for (i in 1:length(plottingConditions)){
  
  if(plottingConditions[i]=="Attended"){dataAmplitudePlot=subset(dataPlot,Attention=="Att")}
  
  if(plottingConditions[i]=="Unattended"){dataAmplitudePlot=subset(dataPlot,Attention=="NotAtt")}  

# Pirate plot

    pirateplot(formula = Amplitude ~ `Reward phase` + `Reward probability`, # dependent~independent variables
             data=dataAmplitudePlot, # data frame
             main=plottingConditions[i], # main title
             ylim=c(0.2,2.2), # y-axis: limits
             ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
             theme=0, # preset theme (0: use your own)
             point.col="black", # points: color
             point.o=.3, # points: opacity (0-1)
             avg.line.col="black", # average line: color
             avg.line.lwd=2, # average line: line width
             avg.line.o=1, # average line: opacity (0-1)
             bean.b.col="black", # bean border, color
             bean.lwd=0.6, # bean border, line width
             bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
             bean.b.o=0.3, # bean border, opacity (0-1)
             bean.f.col="gray", # bean filling, color
             bean.f.o=.1, # bean filling, opacity (0-1)
             cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
             gl.col="gray", # gridlines: color
             gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
             cex.lab=1, # axis labels: size
             cex.axis=1, # axis numbers: size
             cex.names = 1,
             bty="l", # plot box type
             back.col="white") # background, color
}

Statistics

# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))

model.rewardTimesPhasePlusAtt = readRDS("rewardTimesPhasePlusAtt.EEG.alltrials.rds")
model.full = readRDS("full.EEG.alltrials.rds")
compare.threefactors.EEG.waic = readRDS("compare.EEG.waic.alltrials.rds")
bR2.null.EEG = readRDS("bR2.null.EEG.alltrials")
bR2.expphase.EEG = readRDS("bR2.expphase.EEG.alltrials")
bR2.attention.EEG = readRDS("bR2.attention.EEG.alltrials")
bR2.phaseANDattention.EEG = readRDS("bR2.phaseANDattention.EEG.alltrials")
bR2.phaseANDattention_interaction.EEG = readRDS("bR2.phaseANDattention_interaction.EEG.alltrials")
bR2.rewardTimesPhasePlusAtt.EEG = readRDS("bR2.rewardTimesPhasePlusAtt.EEG.alltrials")
bR2.full.EEG = readRDS("bR2.full.EEG.alltrials")

Model comparison with WAIC

print(compare.threefactors.EEG.waic)
## Output of model 'null':
## 
## Computed from 16000 by 1032 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic     11.1 28.1
## p_waic        12.5  0.8
## waic         -22.1 56.2
## 
## Output of model 'expphase':
## 
## Computed from 16000 by 1032 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic     15.7 27.5
## p_waic        35.4  3.1
## waic         -31.3 55.1
## 
## Output of model 'attention':
## 
## Computed from 16000 by 1032 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic    218.0 33.2
## p_waic        59.0  3.9
## waic        -436.0 66.4
## 
## Output of model 'phaseANDattention':
## 
## Computed from 16000 by 1032 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic    232.6 32.5
## p_waic        82.3  6.0
## waic        -465.2 64.9
## 
## Output of model 'phaseANDattention_interaction':
## 
## Computed from 16000 by 1032 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic    230.5 32.6
## p_waic        84.5  6.1
## waic        -461.1 65.2
## 
## Output of model 'rewardTimesPhasePlusAtt':
## 
## Computed from 16000 by 1032 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic    347.9 36.0
## p_waic       129.9  9.3
## waic        -695.9 72.0
## 
## Output of model 'full':
## 
## Computed from 16000 by 1032 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic    344.9 35.9
## p_waic       134.0  9.5
## waic        -689.8 71.8
## 
## Model comparisons:
##                               elpd_diff se_diff
## rewardTimesPhasePlusAtt          0.0       0.0 
## full                            -3.0       2.3 
## phaseANDattention             -115.3      18.0 
## phaseANDattention_interaction -117.4      18.1 
## attention                     -129.9      19.3 
## expphase                      -332.3      25.2 
## null                          -336.9      25.4
#Print the bayes factor for the comparison between the two most complex models
# model.full = readRDS("full.EEG.alltrials.rds")


# bayes_factor(model.full,model.rewardTimesPhasePlusAtt)

Bayesian R squared

Null model

print(bR2.null.EEG)
##      Estimate   Est.Error         Q2.5      Q97.5
## R2 0.01159735 0.009911782 3.619917e-05 0.03559651

Experiment phase

print(bR2.expphase.EEG)
##      Estimate  Est.Error       Q2.5     Q97.5
## R2 0.05029447 0.01570713 0.02113642 0.0824229

Attention model

print(bR2.attention.EEG)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.3696902 0.02058989 0.3286679 0.4091646

Phase and attention model

print(bR2.phaseANDattention.EEG)
##     Estimate  Est.Error      Q2.5    Q97.5
## R2 0.4055565 0.01998488 0.3654231 0.443295

Phase and attention interaction model

print(bR2.phaseANDattention_interaction.EEG)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.4059535 0.02012843 0.3649755 0.4439101

Reward times phase plus attention model

print(bR2.rewardTimesPhasePlusAtt.EEG)
##     Estimate Est.Error      Q2.5     Q97.5
## R2 0.5504777 0.0160925 0.5175255 0.5808094

Full model

print(bR2.full.EEG)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.5521078 0.01611302 0.5192417 0.5827281

Checking the best model

Plotting the chains

# Plot chains
plot(model.rewardTimesPhasePlusAtt, pars = "^b_", ask = FALSE, N=6)

Summary of the best model

# Summary of the best model
print(tidy(model.rewardTimesPhasePlusAtt, par_type = "non-varying"), digits = 1)
##                           term estimate std.error  lower  upper
## 1                    Intercept    1.121      0.02  1.087  1.157
## 2             ConditionLow_Rew    0.001      0.03 -0.052  0.055
## 3                  ExpPhaseAcq    0.025      0.02 -0.006  0.054
## 4                  ExpPhaseExt   -0.007      0.02 -0.045  0.031
## 5              AttentionNotAtt   -0.244      0.03 -0.286 -0.203
## 6 ConditionLow_Rew:ExpPhaseAcq   -0.036      0.02 -0.076  0.005
## 7 ConditionLow_Rew:ExpPhaseExt    0.014      0.02 -0.027  0.054

Plotting the best model

post = posterior_samples(model.rewardTimesPhasePlusAtt, "^b")

# Calculate posteriors for each condition

################################################ Baseline ####

##################### Attended

######### High reward
Baseline_High_Attended = post[["b_Intercept"]]
######### Low reward
Baseline_Low_Attended = post[["b_Intercept"]] + 
  post[["b_ConditionLow_Rew"]] 

##################### Not Attended

######### High reward
Baseline_High_NotAttended = post[["b_Intercept"]] + 
  post[["b_AttentionNotAtt"]]
######### Low reward
Baseline_Low_NotAttended = post[["b_Intercept"]] + 
  post[["b_AttentionNotAtt"]] + 
  post[["b_ConditionLow_Rew"]] 

################################################ Acquistion

##################### Attended

######### High reward
Acquisition_High_Attended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] 
######### Low reward
Acquisition_Low_Attended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ConditionLow_Rew:ExpPhaseAcq"]]

##################### Not Attended

######### High reward
Acquisition_High_NotAttended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] + 
  post[["b_AttentionNotAtt"]] 
  
######### Low reward
Acquisition_Low_NotAttended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] + 
  post[["b_AttentionNotAtt"]] + 
  post[["b_ConditionLow_Rew"]] +
  post[["b_ConditionLow_Rew:ExpPhaseAcq"]] 


################################################ Extinction

##################### Attended

######### High reward
Extinction_High_Attended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] 
######### Low reward
Extinction_Low_Attended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ConditionLow_Rew:ExpPhaseExt"]]

##################### Not Attended

######### High reward
Extinction_High_NotAttended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] + 
  post[["b_AttentionNotAtt"]]
######### Low reward
Extinction_Low_NotAttended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] + 
  post[["b_AttentionNotAtt"]] + 
  post[["b_ConditionLow_Rew"]] +
  post[["b_ConditionLow_Rew:ExpPhaseExt"]] 
# Data from the model

# make a data frame of conditions from the posterior

posterior_conditions = melt(data.frame(Baseline_High_Attended, Baseline_High_NotAttended, Baseline_Low_Attended, Baseline_Low_NotAttended, Acquisition_High_Attended, Acquisition_High_NotAttended, Acquisition_Low_Attended, Acquisition_Low_NotAttended, Extinction_High_Attended, Extinction_High_NotAttended, Extinction_Low_Attended, Extinction_Low_NotAttended))

posterior_conditions =  posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability", "Attention"), "_", extra = "merge")

posterior_conditions$Attention = recode(posterior_conditions$Attention,
                           "Attended" = "Attended",
                           "NotAttended" = "Unattended")

posterior_conditions$Attention  = as.factor(posterior_conditions$Attention )

posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
                                      "High" = "High reward probability",
                                      "Low" = "Low reward probability")

posterior_conditions$`Reward Probability` = as.factor(posterior_conditions$`Reward Probability`)

posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
                                "Acquisition" = "Training",
                                "Baseline" = "Baseline",
                                "Extinction" = "Test")

posterior_conditions$`Reward Phase` = as.factor(posterior_conditions$`Reward Phase`)

names(posterior_conditions)[4] = "Amplitude"




# Make a table with credible intervals
posterior_means = as.data.frame(c("Attended Baseline High Reward", 
                                  "Unattended Baseline High Reward", 
                                  "Attended Baseline Low Reward",
                                  "Unattended Baseline Low Reward",
                                  
                                  "Attended Acquisition High Reward", 
                                  "Unattended Acquisition High Reward", 
                                  "Attended Acquisition Low Reward",
                                  "Unattended Acquisition Low Reward",
                                  
                                  "Attended Extinction High Reward", 
                                  "Unattended Extinction High Reward", 
                                  "Attended Extinction Low Reward",
                                  "Unattended Extinction Low Reward"))

names(posterior_means)[1] = "Condition"

posterior_means$low_95_CI = c(round(hdi(Baseline_High_Attended)[[1]], digits = 2),
                              round(hdi(Baseline_High_NotAttended)[[1]], digits = 2),
                              round(hdi(Baseline_Low_Attended)[[1]], digits = 2),
                              round(hdi(Baseline_Low_NotAttended)[[1]], digits = 2),
                              round(hdi(Acquisition_High_Attended)[[1]], digits = 2),
                              round(hdi(Acquisition_High_NotAttended)[[1]], digits = 2),
                              round(hdi(Acquisition_Low_Attended)[[1]], digits = 2),
                              round(hdi(Acquisition_Low_NotAttended)[[1]], digits = 2),
                              round(hdi(Extinction_High_Attended)[[1]], digits = 2),
                              round(hdi(Extinction_High_NotAttended)[[1]], digits = 2),
                              round(hdi(Extinction_Low_Attended)[[1]], digits = 2),
                              round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2))

posterior_means$high_95_CI = c(round(hdi(Baseline_High_Attended)[[2]], digits = 2),
                              round(hdi(Baseline_High_NotAttended)[[2]], digits = 2),
                              round(hdi(Baseline_Low_Attended)[[2]], digits = 2),
                              round(hdi(Baseline_Low_NotAttended)[[2]], digits = 2),
                              round(hdi(Acquisition_High_Attended)[[2]], digits = 2),
                              round(hdi(Acquisition_High_NotAttended)[[2]], digits = 2),
                              round(hdi(Acquisition_Low_Attended)[[2]], digits = 2),
                              round(hdi(Acquisition_Low_NotAttended)[[2]], digits = 2),
                              round(hdi(Extinction_High_Attended)[[2]], digits = 2),
                              round(hdi(Extinction_High_NotAttended)[[2]], digits = 2),
                              round(hdi(Extinction_Low_Attended)[[2]], digits = 2),
                              round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2))

posterior_means$low_50_CI = c(round(hdi(Baseline_High_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Baseline_High_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Baseline_Low_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Baseline_Low_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Acquisition_High_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Acquisition_High_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Acquisition_Low_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Acquisition_Low_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Extinction_High_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Extinction_High_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Extinction_Low_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Extinction_Low_NotAttended,0.5)[[1]], digits = 2))

posterior_means$high_50_CI = c(round(hdi(Baseline_High_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Baseline_High_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Baseline_Low_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Baseline_Low_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Acquisition_High_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Acquisition_High_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Acquisition_Low_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Acquisition_Low_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Extinction_High_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Extinction_High_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Extinction_Low_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Extinction_Low_NotAttended,0.5)[[2]], digits = 2))


posterior_means =  posterior_means %>% separate(Condition, c("Attention","Reward Phase", "Reward Probability"), " ", extra = "merge")

posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
                                      "High Reward" = "High reward probability",
                                      "Low Reward" = "Low reward probability")

posterior_means$`Reward Probability`= as.factor(posterior_means$`Reward Probability`)

posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
                                "Acquisition" = "Training",
                                "Baseline" = "Baseline",
                                "Extinction" = "Test")

posterior_means$`Reward Phase` = as.factor(posterior_means$`Reward Phase`)

posterior_means$Attention = recode(posterior_means$Attention,
                           "Attended" = "Attended",
                           "Unattended" = "Unattended")

posterior_means$Attention = as.factor(posterior_means$Attention)

# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$Amplitude = c(1,1,1,1,1,1,1,1,1,1,1,1)


# Raw data
# Prepare the dataset
# prepare data for plotting - average over motion

dataPlot = data_all
dataPlot = ddply(dataPlot,.(Attention,ExpPhase,Condition,Subject),plyr::summarize,Amplitude=mean(Amplitude, na.rm = TRUE))


# rename variables
colnames(dataPlot)[colnames(dataPlot)=="ExpPhase"] <- "Reward Phase"
colnames(dataPlot)[colnames(dataPlot)=="Condition"] <- "Reward Probability"

# rename conditions
dataPlot$`Reward Phase` = recode(dataPlot$`Reward Phase`,
                                  "Acq" = "Training",
                                  "Bsln" = "Baseline",
                                  "Ext" = "Test")

dataPlot$`Reward Probability` = recode(dataPlot$`Reward Probability`,
                                        "High_Rew" = "High reward probability",
                                        "Low_Rew" = "Low reward probability")

dataPlot$Attention = recode(dataPlot$Attention,
                                        "Att" = "Attended",
                                        "NotAtt" = "Unattended")


##### Plotting targets and distractors separately  #####
  
# # Attended
#   
#   dataAmplitudePlot=subset(dataPlot,Attention=="Attended")
#   posterior_means_plot=subset(posterior_means,Attention=="Attended")
#   posterior_conditions_plot=subset(posterior_conditions,Attention=="Attended")
#   
#   # tiff("ssvep_attended.tiff", units="in", width=8, height=5, res=800)
#     
#   ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) + 
#     # violin of the raw data
#     geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
#     # individual participants for the raw data
#     # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
#     geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
#     # mean of the model
#     stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
#     # 95 credible intervals of the model
#     geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
#     # 50 credible intervals of the model
#     geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
#     scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
#     facet_grid(.~`Reward Probability`) + 
#     scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
#     theme(
#       axis.line = element_line(size=1, colour = "black"),
#       panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
#       panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
#       panel.border = element_blank(), 
#       panel.background = element_blank(),
#       plot.title = element_text(size = 12,  face = "bold",hjust = 0.5),
#       text=element_text(colour="black", size = 11),
#       axis.text.x=element_text(colour="black", size = 11),
#       axis.text.y=element_text(colour="black", size = 11),
#       axis.title=element_text(size=12,colour = "black")) + 
#     theme(strip.background =element_rect(fill="white")) + 
#     theme(strip.text = element_text(size = 12))
#   
#   # dev.off()
#   
#   
#   
#   # Unattended
#   
#   dataAmplitudePlot=subset(dataPlot,Attention=="Unattended")
#   posterior_means_plot=subset(posterior_means,Attention=="Unattended")
#   posterior_conditions_plot=subset(posterior_conditions,Attention=="Unattended")
#   
#   # tiff("ssvep_unattended.tiff", units="in", width=8, height=5, res=800)
#     
#   ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) + 
#     # violin of the raw data
#     geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
#     # individual participants for the raw data
#     # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
#     geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
#     # mean of the model
#     stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
#     # 95 credible intervals of the model
#     geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
#     # 50 credible intervals of the model
#     geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
#     scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
#     facet_grid(.~`Reward Probability`) + 
#     scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
#     theme(
#       axis.line = element_line(size=1, colour = "black"),
#       panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
#       panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
#       panel.border = element_blank(), 
#       panel.background = element_blank(),
#       plot.title = element_text(size = 12,  face = "bold",hjust = 0.5),
#       text=element_text(colour="black", size = 11),
#       axis.text.x=element_text(colour="black", size = 11),
#       axis.text.y=element_text(colour="black", size = 11),
#       axis.title=element_text(size=12,colour = "black")) + 
#     theme(strip.background =element_rect(fill="white")) + 
#     theme(strip.text = element_text(size = 12))
#   
#   # dev.off()

  
  
##### Plotting high and low separately  #####
  
# High
  
  dataAmplitudePlot=subset(dataPlot,`Reward Probability`=="High reward probability")
  posterior_means_plot=subset(posterior_means,`Reward Probability`=="High reward probability")
  posterior_conditions_plot=subset(posterior_conditions,`Reward Probability`=="High reward probability")
  
  dataPlot$Attention = relevel(dataPlot$Attention,"Attended")
  posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Attended")
  posterior_means$Attention = relevel(posterior_means$Attention,"Attended")
  
  # tiff("ssvep_high.tiff", units="in", width=8, height=5, res=800)
    
f.2.c =   ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=Attention)) + 
    # violin of the raw data
    geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
    # individual participants for the raw data
    # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
    geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
    # mean of the model
    stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=Attention),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
    # 95 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
    # 50 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
    scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
    facet_grid(.~Attention) + 
    scale_y_continuous(name="Amplitude (a.u.)",limits = c(0.6,1.75),breaks=seq(0.6,1.75,0.25)) +
    labs(title="High reward") +
  theme(
    axis.line = element_line(size=1, colour = "black"),
    panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
    panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
    panel.border = element_blank(), 
    panel.background = element_blank(),
    plot.title = element_text(size = 18,  face = "bold",hjust = 0.5),
    text=element_text(colour="black", size = 14),
    axis.text.x=element_text(colour="black", size = 14),
    axis.text.y=element_text(colour="black", size = 14),
    axis.title=element_text(size=16,colour = "black")) + 
  theme(strip.background =element_rect(fill="white")) + 
  theme(strip.text = element_text(size = 16))
f.2.c  

  # dev.off()
  
  
  
# Low
  
  dataAmplitudePlot=subset(dataPlot,`Reward Probability`=="Low reward probability")
  posterior_means_plot=subset(posterior_means,`Reward Probability`=="Low reward probability")
  posterior_conditions_plot=subset(posterior_conditions,`Reward Probability`=="Low reward probability")
  
  dataPlot$Attention = relevel(dataPlot$Attention,"Attended")
  posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Attended")
  posterior_means$Attention = relevel(posterior_means$Attention,"Attended")
  
  # tiff("ssvep_low.tiff", units="in", width=8, height=5, res=800)
    
f.2.d =  ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=Attention)) + 
    # violin of the raw data
    geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
    # individual participants for the raw data
    # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
    geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
    # mean of the model
    stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=Attention),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
    # 95 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
    # 50 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
    scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
    facet_grid(.~Attention) + 
    scale_y_continuous(name="Amplitude (a.u.)",limits = c(0.6,1.75),breaks=seq(0.6,1.75,0.25)) +
    labs(title="Low reward") +
  theme(
    axis.line = element_line(size=1, colour = "black"),
    panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
    panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
    panel.border = element_blank(), 
    panel.background = element_blank(),
    plot.title = element_text(size = 18,  face = "bold",hjust = 0.5),
    text=element_text(colour="black", size = 14),
    axis.text.x=element_text(colour="black", size = 14),
    axis.text.y=element_text(colour="black", size = 14),
    axis.title=element_text(size=16,colour = "black")) + 
  theme(strip.background =element_rect(fill="white")) + 
  theme(strip.text = element_text(size = 16))
f.2.d  

  # dev.off()
####### Plotting everything together #####
  
# # Order for plotting
# # dataPlot$Attention = factor(dataPlot$Attention, levels = c("Target","Distractor")) 
# # posterior_conditions$Attention = factor(posterior_conditions$Attention, levels = c("Target","Distractor")) 
# # posterior_means$Attention = factor(posterior_means$Attention, levels = c("Target","Distractor")) 
# # 
# # dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward Phase`,dataPlot$`Reward Probability`),]
# 
# dataPlot$Attention = relevel(dataPlot$Attention,"Attended")
# posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Attended")
# posterior_means$Attention = relevel(posterior_means$Attention,"Attended")
# 
# 
#   
#   
#     # tiff("ssvep_both.tiff", units="in", width=16, height=10, res=800)
# 
#     ggplot(data=dataPlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) + 
#     # violin of the raw data
#     geom_violin(data=dataPlot, color ="#636363",trim = TRUE,alpha="0.2") +
#     # individual participants for the raw data
#     # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
#     geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
#     # mean of the model
#     stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
#     # 95 credible intervals of the model
#     geom_linerange(data = posterior_means,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
#     # 50 credible intervals of the model
#     geom_linerange(data = posterior_means,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
#     scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
#     facet_grid(`Reward Probability`~Attention) + 
#     scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(0.5,2,0.5)) +
#     theme(
#       axis.line = element_line(size=2, colour = "black"),
#       panel.grid.major = element_line(colour = "grey",size = 0.2,linetype = "dashed"),
#       panel.grid.minor = element_line(colour = "grey",size = 0.2,linetype = "dashed"),
#       panel.border = element_blank(), 
#       panel.background = element_blank(),
#       plot.title = element_text(size = 24,  face = "bold",hjust = 0.5),
#       text=element_text(colour="black", size = 22),
#       axis.text.x=element_text(colour="black", size = 22),
#       axis.text.y=element_text(colour="black", size = 22),
#       axis.title=element_text(size=24,colour = "black")) + 
#     theme(strip.background =element_rect(fill="white")) + 
#     theme(strip.text = element_text(size = 24))
# # dev.off()

Table of means across conditions

# Make a table with conditions
posterior_means = as.data.frame(c("Attended Baseline High Reward", 
                                  "Attended Baseline Low Reward", 
                                  "Attended Acquisition High Reward", 
                                  "Attended Acquisition Low Reward", 
                                  "Attended Extinction High Reward", 
                                  "Attended Extinction Low Reward",
                                  "Unattended Baseline High Reward", 
                                  "Unattended Baseline Low Reward", 
                                  "Unattended Acquisition High Reward", 
                                  "Unattended Acquisition Low Reward", 
                                  "Unattended Extinction High Reward", 
                                  "Unattended Extinction Low Reward"))
names(posterior_means)[1] = "Condition"

posterior_means$Mean = c(paste(round(mean(Baseline_High_Attended), digits = 2), " [", round(hdi(Baseline_High_Attended)[[1]], digits = 2), " ", round(hdi(Baseline_High_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Baseline_Low_Attended), digits = 2), " [", round(hdi(Baseline_Low_Attended)[[1]], digits = 2), " ", round(hdi(Baseline_Low_Attended)[[2]], digits = 2), "]"),
                        
                        paste(round(mean(Acquisition_High_Attended), digits = 2), " [", round(hdi(Acquisition_High_Attended)[[1]], digits = 2), " ", round(hdi(Acquisition_High_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Acquisition_Low_Attended), digits = 2), " [", round(hdi(Acquisition_Low_Attended)[[1]], digits = 2), " ", round(hdi(Acquisition_Low_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_High_Attended), digits = 2), " [", round(hdi(Extinction_High_Attended)[[1]], digits = 2), " ", round(hdi(Extinction_High_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_Low_Attended), digits = 2), " [", round(hdi(Extinction_Low_Attended)[[1]], digits = 2), " ", round(hdi(Extinction_Low_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Baseline_High_NotAttended), digits = 2), " [", round(hdi(Baseline_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Baseline_High_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Baseline_Low_NotAttended), digits = 2), " [", round(hdi(Baseline_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Baseline_Low_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Acquisition_High_NotAttended), digits = 2), " [", round(hdi(Acquisition_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Acquisition_High_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Acquisition_Low_NotAttended), digits = 2), " [", round(hdi(Acquisition_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Acquisition_Low_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_High_NotAttended), digits = 2), " [", round(hdi(Extinction_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_High_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_Low_NotAttended), digits = 2), " [", round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2), "]"))

names(posterior_means)[2] = "Mean [HDI]"

posterior_means =  posterior_means %>% separate(Condition, c("Attention", "Reward Phase", "Reward Probability"), " ", extra = "merge")

kable(posterior_means, caption = "Means per condition")
Means per condition
Attention Reward Phase Reward Probability Mean [HDI]
Attended Baseline High Reward 1.12 [ 1.08 1.16 ]
Attended Baseline Low Reward 1.12 [ 1.08 1.17 ]
Attended Acquisition High Reward 1.15 [ 1.1 1.19 ]
Attended Acquisition Low Reward 1.11 [ 1.06 1.16 ]
Attended Extinction High Reward 1.11 [ 1.06 1.16 ]
Attended Extinction Low Reward 1.13 [ 1.07 1.19 ]
Unattended Baseline High Reward 0.88 [ 0.83 0.92 ]
Unattended Baseline Low Reward 0.88 [ 0.84 0.92 ]
Unattended Acquisition High Reward 0.9 [ 0.85 0.95 ]
Unattended Acquisition Low Reward 0.87 [ 0.83 0.91 ]
Unattended Extinction High Reward 0.87 [ 0.82 0.93 ]
Unattended Extinction Low Reward 0.89 [ 0.83 0.94 ]
### Inference about the best model

Attended vs. unattended

Check the difference between attended and not attended in baseline high rewarded

Diff_Att_NotAtt_Bsln_High = Baseline_High_Attended - Baseline_High_NotAttended
plotPost(Diff_Att_NotAtt_Bsln_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

length(Diff_Att_NotAtt_Bsln_High)
## [1] 16000

Check the difference between attended and not attended in baseline low rewarded

Diff_Att_NotAtt_Bsln_Low = Baseline_Low_Attended - Baseline_Low_NotAttended
plotPost(Diff_Att_NotAtt_Bsln_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

length(Diff_Att_NotAtt_Bsln_Low)
## [1] 16000

Check the difference between attended and not attended in acquisition high rewarded

Diff_Att_NotAtt_Acq_High = Acquisition_High_Attended - Acquisition_High_NotAttended
plotPost(Diff_Att_NotAtt_Acq_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

length(Diff_Att_NotAtt_Acq_High)
## [1] 16000

Check the difference between attended and not attended in acquisition low rewarded

Diff_Att_NotAtt_Acq_Low = Acquisition_Low_Attended - Acquisition_Low_NotAttended
plotPost(Diff_Att_NotAtt_Acq_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

length(Diff_Att_NotAtt_Acq_Low)
## [1] 16000

Check the difference between attended and not attended in extinction high rewarded

Diff_Att_NotAtt_Ext_High = Extinction_High_Attended - Extinction_High_NotAttended
plotPost(Diff_Att_NotAtt_Ext_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

length(Diff_Att_NotAtt_Ext_High)
## [1] 16000

Check the difference between attended and not attended in extinction low rewarded

Diff_Att_NotAtt_Ext_Low = Extinction_Low_Attended - Extinction_Low_NotAttended
plotPost(Diff_Att_NotAtt_Ext_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Comparison between phases

Check the difference between baseline and acquisition in high reward attended

Diff_Bsln_Acq_High_Att = Acquisition_High_Attended - Baseline_High_Attended
plotPost(Diff_Bsln_Acq_High_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between baseline and acquisition in low reward attended

Diff_Bsln_Acq_Low_Att = Baseline_Low_Attended - Acquisition_Low_Attended
plotPost(Diff_Bsln_Acq_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between baseline and acquisition in high reward not attended

Diff_Bsln_Acq_High_NotAtt =  Acquisition_High_NotAttended - Baseline_High_NotAttended
plotPost(Diff_Bsln_Acq_High_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between baseline and acquisition in low reward not attended

Diff_Bsln_Acq_Low_NotAtt = Acquisition_Low_NotAttended - Baseline_Low_NotAttended  
plotPost(Diff_Bsln_Acq_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between acquisition and extinction in high reward attended

Diff_Acq_Ext_High_Att = Acquisition_High_Attended - Extinction_High_Attended
plotPost(Diff_Acq_Ext_High_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between acquisition and extinction in low reward attended

Diff_Acq_Ext_Low_Att = Extinction_Low_Attended - Acquisition_Low_Attended 
plotPost(Diff_Acq_Ext_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between acquisition and extinction in high reward not attended

Diff_Acq_Ext_High_NotAtt = Extinction_High_NotAttended - Acquisition_High_NotAttended 
plotPost(Diff_Acq_Ext_High_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between acquisition and extinction in low reward not attended

Diff_Acq_Ext_Low_NotAtt = Extinction_Low_NotAttended - Acquisition_Low_NotAttended 
plotPost(Diff_Acq_Ext_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Baseline difference

Check the difference between high and low reward in baseline attended

Diff_Bsln_High_Low_Att = Baseline_High_Attended - Baseline_Low_Attended
plotPost(Diff_Bsln_High_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between high and low reward in baseline not attended

Diff_Bsln_High_Low_NotAtt = Baseline_High_NotAttended - Baseline_Low_NotAttended
plotPost(Diff_Bsln_High_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between baseline and acquisition in high reward unattended vs. attended

Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended = Diff_Bsln_High_Low_NotAtt - Diff_Bsln_High_Low_Att 
plotPost(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

paste("Mean = ",round(mean(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended), digits = 2), " [", round(hdi(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended)[[1]], digits = 2), " ", round(hdi(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended)[[2]], digits = 2), "]")
## [1] "Mean =  0  [ 0   0 ]"

Only movement trials

# Take only the movement trials
data = subset(data_all, Movement=="Mov")

Means - raw data

summary = ddply(data,.(Attention,ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(Amplitude, na.rm = TRUE), digits = 2), " [", round(hdi(Amplitude)[[1]], digits = 2), " ", round(hdi(Amplitude)[[2]], digits = 2), "]")))

names(summary) = c("Attention", "Reward phase", "Reward probability", "Amplitude")

summary$Attention = recode(summary$Attention,
                           "Att" = "Attended",
                           "NotAtt" = "Unattended")

summary$`Reward phase` = recode(summary$`Reward phase`,
                           "Acq" = "Acquisition",
                           "Bsln" = "Baseline",
                           "Ext" = "Extinction")

summary$`Reward probability` = recode(summary$`Reward probability`,
                           "High_Rew" = "High",
                           "Low_Rew" = "Low")

summary = as.data.frame(summary)

summary$`Reward phase` = factor(summary$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
summary = summary[order(summary$Attention,summary$`Reward phase`,summary$`Reward probability`),]
row.names(summary) = NULL

kable(summary, caption = "Amplitudes per condition")
Amplitudes per condition
Attention Reward phase Reward probability Amplitude
Attended Baseline High 1.1 [ 0.77 1.39 ]
Attended Baseline Low 1.12 [ 0.86 1.52 ]
Attended Acquisition High 1.16 [ 0.83 1.6 ]
Attended Acquisition Low 1.09 [ 0.87 1.48 ]
Attended Extinction High 1.13 [ 0.61 1.53 ]
Attended Extinction Low 1.09 [ 0.59 1.59 ]
Unattended Baseline High 0.87 [ 0.5 1.13 ]
Unattended Baseline Low 0.87 [ 0.66 1.07 ]
Unattended Acquisition High 0.88 [ 0.59 1.09 ]
Unattended Acquisition Low 0.87 [ 0.69 1.13 ]
Unattended Extinction High 0.85 [ 0.55 1.1 ]
Unattended Extinction Low 0.89 [ 0.44 1.24 ]

Plots - raw data

# Plot amplitude across experiment phases------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# prepare data for plotting
dataPlot = data

# rename variables
colnames(dataPlot)[colnames(dataPlot)=="ExpPhase"] <- "Reward phase"
colnames(dataPlot)[colnames(dataPlot)=="Condition"] <- "Reward probability"

# rename conditions
dataPlot$`Reward phase` = recode(dataPlot$`Reward phase`,
                                  "Acq" = "Acquisition",
                                  "Bsln" = "Baseline",
                                  "Ext" = "Extinction")

dataPlot$`Reward probability` = recode(dataPlot$`Reward probability`,
                                        "High_Rew" = "High",
                                        "Low_Rew" = "Low")

#order
dataPlot$`Reward phase` = factor(dataPlot$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward phase`,dataPlot$`Reward probability`),]


plottingConditions = c("Attended","Unattended" )
for (i in 1:length(plottingConditions)){
  
  if(plottingConditions[i]=="Attended"){dataAmplitudePlot=subset(dataPlot,Attention=="Att")}
  
  if(plottingConditions[i]=="Unattended"){dataAmplitudePlot=subset(dataPlot,Attention=="NotAtt")}  

# Pirate plot

    pirateplot(formula = Amplitude ~ `Reward phase` + `Reward probability`, # dependent~independent variables
             data=dataAmplitudePlot, # data frame
             main=plottingConditions[i], # main title
             ylim=c(0.2,2.2), # y-axis: limits
             ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
             theme=0, # preset theme (0: use your own)
             point.col="black", # points: color
             point.o=.3, # points: opacity (0-1)
             avg.line.col="black", # average line: color
             avg.line.lwd=2, # average line: line width
             avg.line.o=1, # average line: opacity (0-1)
             bean.b.col="black", # bean border, color
             bean.lwd=0.6, # bean border, line width
             bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
             bean.b.o=0.3, # bean border, opacity (0-1)
             bean.f.col="gray", # bean filling, color
             bean.f.o=.1, # bean filling, opacity (0-1)
             cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
             gl.col="gray", # gridlines: color
             gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
             cex.lab=1, # axis labels: size
             cex.axis=1, # axis numbers: size
             cex.names = 1,
             bty="l", # plot box type
             back.col="white") # background, color
}

Statistics

# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))

model.rewardTimesPhasePlusAtt = readRDS("rewardTimesPhasePlusAtt.EEG.movementtrials.rds")
compare.threefactors.EEG.waic = readRDS("compare.EEG.waic.movementtrials.rds")
bR2.null.EEG = readRDS("bR2.null.EEG.movementtrials")
bR2.expphase.EEG = readRDS("bR2.expphase.EEG.movementtrials")
bR2.attention.EEG = readRDS("bR2.attention.EEG.movementtrials")
bR2.phaseANDattention.EEG = readRDS("bR2.phaseANDattention.EEG.movementtrials")
bR2.phaseANDattention_interaction.EEG = readRDS("bR2.phaseANDattention_interaction.EEG.movementtrials")
bR2.rewardTimesPhasePlusAtt.EEG = readRDS("bR2.rewardTimesPhasePlusAtt.EEG.movementtrials")
bR2.full.EEG = readRDS("bR2.full.EEG.movementtrials")

Model comparison with WAIC

print(compare.threefactors.EEG.waic)
## Output of model 'null':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic     18.5 19.9
## p_waic         6.2  0.7
## waic         -37.0 39.8
## 
## Output of model 'expphase':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic     14.8 19.4
## p_waic        23.6  3.3
## waic         -29.6 38.8
## 
## Output of model 'attention':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic    107.4 25.5
## p_waic        36.2  3.8
## waic        -214.9 51.1
## 
## Output of model 'phaseANDattention':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic    107.9 23.6
## p_waic        61.4  6.7
## waic        -215.8 47.2
## 
## Output of model 'phaseANDattention_interaction':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic    105.2 23.7
## p_waic        62.8  6.8
## waic        -210.4 47.4
## 
## Output of model 'rewardTimesPhasePlusAtt':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic    171.6 24.4
## p_waic       113.6 10.4
## waic        -343.2 48.9
## 
## Output of model 'full':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic    169.2 24.1
## p_waic       117.1 10.4
## waic        -338.3 48.2
## 
## Model comparisons:
##                               elpd_diff se_diff
## rewardTimesPhasePlusAtt          0.0       0.0 
## full                            -2.5       2.8 
## phaseANDattention              -63.7      12.6 
## attention                      -64.2      12.8 
## phaseANDattention_interaction  -66.4      12.6 
## null                          -153.1      15.8 
## expphase                      -156.8      16.2

Bayesian R squared

Null model

print(bR2.null.EEG)
##       Estimate   Est.Error         Q2.5      Q97.5
## R2 0.007491801 0.009307113 9.034355e-06 0.03367014

Experiment phase

print(bR2.expphase.EEG)
##      Estimate Est.Error        Q2.5      Q97.5
## R2 0.04206872  0.022931 0.007001546 0.09329268

Attention model

print(bR2.attention.EEG)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.3399462 0.03301465 0.2726915 0.4016646

Phase and attention model

print(bR2.phaseANDattention.EEG)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.3836096 0.03609204 0.3084417 0.4495137

Phase and attention interaction model

print(bR2.phaseANDattention_interaction.EEG)
##     Estimate Est.Error      Q2.5     Q97.5
## R2 0.3825971 0.0360869 0.3080549 0.4500715

Reward times phase plus attention model

print(bR2.rewardTimesPhasePlusAtt.EEG)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.5781837 0.02530734 0.5245946 0.6234436

Full model

print(bR2.full.EEG)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.5816829 0.02495645 0.5295854 0.6266049

Checking the best model

Plotting the chains

# Plot chains
plot(model.rewardTimesPhasePlusAtt, pars = "^b_", ask = FALSE, N=6)

Summary of the best model

# Summary of the best model
print(tidy(model.rewardTimesPhasePlusAtt, par_type = "non-varying"), digits = 1)
##                           term estimate std.error  lower  upper
## 1                    Intercept    1.106      0.03  1.063  1.149
## 2             ConditionLow_Rew    0.012      0.04 -0.047  0.070
## 3                  ExpPhaseAcq    0.038      0.02 -0.002  0.077
## 4                  ExpPhaseExt    0.009      0.03 -0.041  0.058
## 5              AttentionNotAtt   -0.245      0.02 -0.283 -0.207
## 6 ConditionLow_Rew:ExpPhaseAcq   -0.057      0.03 -0.111 -0.002
## 7 ConditionLow_Rew:ExpPhaseExt   -0.012      0.03 -0.065  0.043

Plotting the best model

post = posterior_samples(model.rewardTimesPhasePlusAtt, "^b")

# Calculate posteriors for each condition

################################################ Baseline ####

##################### Attended

######### High reward
Baseline_High_Attended = post[["b_Intercept"]]
######### Low reward
Baseline_Low_Attended = post[["b_Intercept"]] + 
  post[["b_ConditionLow_Rew"]] 

##################### Not Attended

######### High reward
Baseline_High_NotAttended = post[["b_Intercept"]] + 
  post[["b_AttentionNotAtt"]]
######### Low reward
Baseline_Low_NotAttended = post[["b_Intercept"]] + 
  post[["b_AttentionNotAtt"]] + 
  post[["b_ConditionLow_Rew"]] 

################################################ Acquistion

##################### Attended

######### High reward
Acquisition_High_Attended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] 
######### Low reward
Acquisition_Low_Attended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ConditionLow_Rew:ExpPhaseAcq"]]

##################### Not Attended

######### High reward
Acquisition_High_NotAttended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] + 
  post[["b_AttentionNotAtt"]] 

######### Low reward
Acquisition_Low_NotAttended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] + 
  post[["b_AttentionNotAtt"]] + 
  post[["b_ConditionLow_Rew"]] +
  post[["b_ConditionLow_Rew:ExpPhaseAcq"]] 


################################################ Extinction

##################### Attended

######### High reward
Extinction_High_Attended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] 
######### Low reward
Extinction_Low_Attended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ConditionLow_Rew:ExpPhaseExt"]]

##################### Not Attended

######### High reward
Extinction_High_NotAttended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] + 
  post[["b_AttentionNotAtt"]]
######### Low reward
Extinction_Low_NotAttended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] + 
  post[["b_AttentionNotAtt"]] + 
  post[["b_ConditionLow_Rew"]] +
  post[["b_ConditionLow_Rew:ExpPhaseExt"]] 
# make a data frame

posterior_conditions = melt(data.frame(Baseline_High_Attended, Baseline_High_NotAttended, Baseline_Low_Attended, Baseline_Low_NotAttended, Acquisition_High_Attended, Acquisition_High_NotAttended, Acquisition_Low_Attended, Acquisition_Low_NotAttended, Extinction_High_Attended, Extinction_High_NotAttended, Extinction_Low_Attended, Extinction_Low_NotAttended))

posterior_conditions =  posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability", "Attention"), "_", extra = "merge")

posterior_conditions$Attention = recode(posterior_conditions$Attention,
                           "Attended" = "Attended",
                           "NotAttended" = "Unattended")

names(posterior_conditions)[4] = "Amplitude"


#order
#dataPlot$`Reward phase` = factor(dataPlot$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
#dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward phase`,dataPlot$`Reward probability`),]


plottingConditions = c("Attended","Unattended" )
for (i in 1:length(plottingConditions)){
  
  if(plottingConditions[i]=="Attended"){dataAmplitudePlot=subset(posterior_conditions,Attention=="Attended")}
  
  if(plottingConditions[i]=="Unattended"){dataAmplitudePlot=subset(posterior_conditions,Attention=="Unattended")}  

# Pirate plot

    pirateplot(formula = Amplitude ~ `Reward Phase` + `Reward Probability`, # dependent~independent variables
             data=dataAmplitudePlot, # data frame
             main=plottingConditions[i], # main title
             ylim=c(0.7,1.3), # y-axis: limits
             ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
             theme=0, # preset theme (0: use your own)
             avg.line.col="black", # average line: color
             avg.line.lwd=2, # average line: line width
             avg.line.o=1, # average line: opacity (0-1)
             bean.b.col="black", # bean border, color
             bean.lwd=0.6, # bean border, line width
             bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
             bean.b.o=0.3, # bean border, opacity (0-1)
             bean.f.col="gray", # bean filling, color
             bean.f.o=.1, # bean filling, opacity (0-1)
             cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
             gl.col="gray", # gridlines: color
             gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
             cex.lab=1, # axis labels: size
             cex.axis=1, # axis numbers: size
             cex.names = 1,
             sortx = "sequential",
             bty="l", # plot box type
             back.col="white") # background, color
}

plottingConditions = c("High","Low" )
for (i in 1:length(plottingConditions)){
  
  if(plottingConditions[i]=="High"){dataAmplitudePlot=subset(posterior_conditions,`Reward Probability`=="High")}
  
  if(plottingConditions[i]=="Low"){dataAmplitudePlot=subset(posterior_conditions,`Reward Probability`=="Low")}  

# Pirate plot

    pirateplot(formula = Amplitude ~ `Attention` + `Reward Phase`, # dependent~independent variables
             data=dataAmplitudePlot, # data frame
             main=plottingConditions[i], # main title
             ylim=c(0.7,1.3), # y-axis: limits
             ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
             theme=0, # preset theme (0: use your own)
             avg.line.col="black", # average line: color
             avg.line.lwd=2, # average line: line width
             avg.line.o=1, # average line: opacity (0-1)
             bean.b.col="black", # bean border, color
             bean.lwd=0.6, # bean border, line width
             bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
             bean.b.o=0.3, # bean border, opacity (0-1)
             bean.f.col="gray", # bean filling, color
             bean.f.o=.1, # bean filling, opacity (0-1)
             cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
             gl.col="gray", # gridlines: color
             gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
             cex.lab=1, # axis labels: size
             cex.axis=1, # axis numbers: size
             cex.names = 1,
             sortx = "sequential",
             bty="l", # plot box type
             back.col="white") # background, color
}

# Data from the model

# make a data frame of conditions from the posterior

posterior_conditions = melt(data.frame(Baseline_High_Attended, Baseline_High_NotAttended, Baseline_Low_Attended, Baseline_Low_NotAttended, Acquisition_High_Attended, Acquisition_High_NotAttended, Acquisition_Low_Attended, Acquisition_Low_NotAttended, Extinction_High_Attended, Extinction_High_NotAttended, Extinction_Low_Attended, Extinction_Low_NotAttended))

posterior_conditions =  posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability", "Attention"), "_", extra = "merge")

posterior_conditions$Attention = recode(posterior_conditions$Attention,
                           "Attended" = "Target",
                           "NotAttended" = "Distractor")

posterior_conditions$Attention  = as.factor(posterior_conditions$Attention )

posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
                                      "High" = "High reward probability",
                                      "Low" = "Low reward probability")

posterior_conditions$`Reward Probability` = as.factor(posterior_conditions$`Reward Probability`)

posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
                                "Acquisition" = "Training",
                                "Baseline" = "Baseline",
                                "Extinction" = "Test")

posterior_conditions$`Reward Phase` = as.factor(posterior_conditions$`Reward Phase`)

names(posterior_conditions)[4] = "Amplitude"




# Make a table with credible intervals
posterior_means = as.data.frame(c("Attended Baseline High Reward", 
                                  "Unattended Baseline High Reward", 
                                  "Attended Baseline Low Reward",
                                  "Unattended Baseline Low Reward",
                                  
                                  "Attended Acquisition High Reward", 
                                  "Unattended Acquisition High Reward", 
                                  "Attended Acquisition Low Reward",
                                  "Unattended Acquisition Low Reward",
                                  
                                  "Attended Extinction High Reward", 
                                  "Unattended Extinction High Reward", 
                                  "Attended Extinction Low Reward",
                                  "Unattended Extinction Low Reward"))

names(posterior_means)[1] = "Condition"

posterior_means$low_95_CI = c(round(hdi(Baseline_High_Attended)[[1]], digits = 2),
                              round(hdi(Baseline_High_NotAttended)[[1]], digits = 2),
                              round(hdi(Baseline_Low_Attended)[[1]], digits = 2),
                              round(hdi(Baseline_Low_NotAttended)[[1]], digits = 2),
                              round(hdi(Acquisition_High_Attended)[[1]], digits = 2),
                              round(hdi(Acquisition_High_NotAttended)[[1]], digits = 2),
                              round(hdi(Acquisition_Low_Attended)[[1]], digits = 2),
                              round(hdi(Acquisition_Low_NotAttended)[[1]], digits = 2),
                              round(hdi(Extinction_High_Attended)[[1]], digits = 2),
                              round(hdi(Extinction_High_NotAttended)[[1]], digits = 2),
                              round(hdi(Extinction_Low_Attended)[[1]], digits = 2),
                              round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2))

posterior_means$high_95_CI = c(round(hdi(Baseline_High_Attended)[[2]], digits = 2),
                              round(hdi(Baseline_High_NotAttended)[[2]], digits = 2),
                              round(hdi(Baseline_Low_Attended)[[2]], digits = 2),
                              round(hdi(Baseline_Low_NotAttended)[[2]], digits = 2),
                              round(hdi(Acquisition_High_Attended)[[2]], digits = 2),
                              round(hdi(Acquisition_High_NotAttended)[[2]], digits = 2),
                              round(hdi(Acquisition_Low_Attended)[[2]], digits = 2),
                              round(hdi(Acquisition_Low_NotAttended)[[2]], digits = 2),
                              round(hdi(Extinction_High_Attended)[[2]], digits = 2),
                              round(hdi(Extinction_High_NotAttended)[[2]], digits = 2),
                              round(hdi(Extinction_Low_Attended)[[2]], digits = 2),
                              round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2))

posterior_means$low_50_CI = c(round(hdi(Baseline_High_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Baseline_High_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Baseline_Low_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Baseline_Low_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Acquisition_High_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Acquisition_High_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Acquisition_Low_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Acquisition_Low_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Extinction_High_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Extinction_High_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Extinction_Low_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Extinction_Low_NotAttended,0.5)[[1]], digits = 2))

posterior_means$high_50_CI = c(round(hdi(Baseline_High_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Baseline_High_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Baseline_Low_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Baseline_Low_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Acquisition_High_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Acquisition_High_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Acquisition_Low_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Acquisition_Low_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Extinction_High_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Extinction_High_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Extinction_Low_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Extinction_Low_NotAttended,0.5)[[2]], digits = 2))


posterior_means =  posterior_means %>% separate(Condition, c("Attention","Reward Phase", "Reward Probability"), " ", extra = "merge")

posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
                                      "High Reward" = "High reward probability",
                                      "Low Reward" = "Low reward probability")

posterior_means$`Reward Probability`= as.factor(posterior_means$`Reward Probability`)

posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
                                "Acquisition" = "Training",
                                "Baseline" = "Baseline",
                                "Extinction" = "Test")

posterior_means$`Reward Phase` = as.factor(posterior_means$`Reward Phase`)

posterior_means$Attention = recode(posterior_means$Attention,
                           "Attended" = "Target",
                           "Unattended" = "Distractor")

posterior_means$Attention = as.factor(posterior_means$Attention)

# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$Amplitude = c(1,1,1,1,1,1,1,1,1,1,1,1)


# Raw data
# Prepare the dataset
# prepare data for plotting - average over motion

dataPlot = data_all
dataPlot = ddply(dataPlot,.(Attention,ExpPhase,Condition,Subject),plyr::summarize,Amplitude=mean(Amplitude, na.rm = TRUE))


# rename variables
colnames(dataPlot)[colnames(dataPlot)=="ExpPhase"] <- "Reward Phase"
colnames(dataPlot)[colnames(dataPlot)=="Condition"] <- "Reward Probability"

# rename conditions
dataPlot$`Reward Phase` = recode(dataPlot$`Reward Phase`,
                                  "Acq" = "Training",
                                  "Bsln" = "Baseline",
                                  "Ext" = "Test")

dataPlot$`Reward Probability` = recode(dataPlot$`Reward Probability`,
                                        "High_Rew" = "High reward probability",
                                        "Low_Rew" = "Low reward probability")

dataPlot$Attention = recode(dataPlot$Attention,
                                        "Att" = "Target",
                                        "NotAtt" = "Distractor")


##### Plotting targets and distractors separately  #####
  
# Attended
  
  dataAmplitudePlot=subset(dataPlot,Attention=="Target")
  posterior_means_plot=subset(posterior_means,Attention=="Target")
  posterior_conditions_plot=subset(posterior_conditions,Attention=="Target")
  
  # tiff("ssvep_attended.tiff", units="in", width=8, height=5, res=800)
    
  ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) + 
    # violin of the raw data
    geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
    # individual participants for the raw data
    # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
    geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
    # mean of the model
    stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
    # 95 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
    # 50 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
    scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
    facet_grid(.~`Reward Probability`) + 
    scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
    theme(
      axis.line = element_line(size=1, colour = "black"),
      panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.border = element_blank(), 
      panel.background = element_blank(),
      plot.title = element_text(size = 12,  face = "bold",hjust = 0.5),
      text=element_text(colour="black", size = 11),
      axis.text.x=element_text(colour="black", size = 11),
      axis.text.y=element_text(colour="black", size = 11),
      axis.title=element_text(size=12,colour = "black")) + 
    theme(strip.background =element_rect(fill="white")) + 
    theme(strip.text = element_text(size = 12))

  # dev.off()
  
  
  
  # Unattended
  
  dataAmplitudePlot=subset(dataPlot,Attention=="Distractor")
  posterior_means_plot=subset(posterior_means,Attention=="Distractor")
  posterior_conditions_plot=subset(posterior_conditions,Attention=="Distractor")
  
  # tiff("ssvep_unattended.tiff", units="in", width=8, height=5, res=800)
    
  ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) + 
    # violin of the raw data
    geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
    # individual participants for the raw data
    # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
    geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
    # mean of the model
    stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
    # 95 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
    # 50 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
    scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
    facet_grid(.~`Reward Probability`) + 
    scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
    theme(
      axis.line = element_line(size=1, colour = "black"),
      panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.border = element_blank(), 
      panel.background = element_blank(),
      plot.title = element_text(size = 12,  face = "bold",hjust = 0.5),
      text=element_text(colour="black", size = 11),
      axis.text.x=element_text(colour="black", size = 11),
      axis.text.y=element_text(colour="black", size = 11),
      axis.title=element_text(size=12,colour = "black")) + 
    theme(strip.background =element_rect(fill="white")) + 
    theme(strip.text = element_text(size = 12))

  # dev.off()

  
  
##### Plotting high and low separately  #####
  
# High
  
  dataAmplitudePlot=subset(dataPlot,`Reward Probability`=="High reward probability")
  posterior_means_plot=subset(posterior_means,`Reward Probability`=="High reward probability")
  posterior_conditions_plot=subset(posterior_conditions,`Reward Probability`=="High reward probability")
  
  dataPlot$Attention = relevel(dataPlot$Attention,"Target")
  posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
  posterior_means$Attention = relevel(posterior_means$Attention,"Target")
  
  # tiff("ssvep_high.tiff", units="in", width=8, height=5, res=800)
    
  ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=Attention)) + 
    # violin of the raw data
    geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
    # individual participants for the raw data
    # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
    geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
    # mean of the model
    stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=Attention),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
    # 95 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
    # 50 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
    scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
    facet_grid(.~Attention) + 
    scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
    labs(title="High reward") +
    theme(
      axis.line = element_line(size=1, colour = "black"),
      panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.border = element_blank(), 
      panel.background = element_blank(),
      plot.title = element_text(size = 12,  face = "bold",hjust = 0.5),
      text=element_text(colour="black", size = 11),
      axis.text.x=element_text(colour="black", size = 11),
      axis.text.y=element_text(colour="black", size = 11),
      axis.title=element_text(size=12,colour = "black")) + 
    theme(strip.background =element_rect(fill="white")) + 
    theme(strip.text = element_text(size = 12))

  # dev.off()
  
  
  
# Low
  
  dataAmplitudePlot=subset(dataPlot,`Reward Probability`=="Low reward probability")
  posterior_means_plot=subset(posterior_means,`Reward Probability`=="Low reward probability")
  posterior_conditions_plot=subset(posterior_conditions,`Reward Probability`=="Low reward probability")
  
  dataPlot$Attention = relevel(dataPlot$Attention,"Target")
  posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
  posterior_means$Attention = relevel(posterior_means$Attention,"Target")
  
  # tiff("ssvep_low.tiff", units="in", width=8, height=5, res=800)
    
  ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=Attention)) + 
    # violin of the raw data
    geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
    # individual participants for the raw data
    # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
    geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
    # mean of the model
    stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=Attention),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
    # 95 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
    # 50 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
    scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
    facet_grid(.~Attention) + 
    scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
    labs(title="Low reward") +
    theme(
      axis.line = element_line(size=1, colour = "black"),
      panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.border = element_blank(), 
      panel.background = element_blank(),
      plot.title = element_text(size = 12,  face = "bold",hjust = 0.5),
      text=element_text(colour="black", size = 11),
      axis.text.x=element_text(colour="black", size = 11),
      axis.text.y=element_text(colour="black", size = 11),
      axis.title=element_text(size=12,colour = "black")) + 
    theme(strip.background =element_rect(fill="white")) + 
    theme(strip.text = element_text(size = 12))

  # dev.off()
####### Plotting everything together #####
  
# Order for plotting
# dataPlot$Attention = factor(dataPlot$Attention, levels = c("Target","Distractor")) 
# posterior_conditions$Attention = factor(posterior_conditions$Attention, levels = c("Target","Distractor")) 
# posterior_means$Attention = factor(posterior_means$Attention, levels = c("Target","Distractor")) 
# 
# dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward Phase`,dataPlot$`Reward Probability`),]

dataPlot$Attention = relevel(dataPlot$Attention,"Target")
posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
posterior_means$Attention = relevel(posterior_means$Attention,"Target")


  
  
    # tiff("ssvep_both.tiff", units="in", width=16, height=10, res=800)

    ggplot(data=dataPlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) + 
    # violin of the raw data
    geom_violin(data=dataPlot, color ="#636363",trim = TRUE,alpha="0.2") +
    # individual participants for the raw data
    # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
    geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
    # mean of the model
    stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
    # 95 credible intervals of the model
    geom_linerange(data = posterior_means,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
    # 50 credible intervals of the model
    geom_linerange(data = posterior_means,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
    scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
    facet_grid(`Reward Probability`~Attention) + 
    scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(0.5,2,0.5)) +
    theme(
      axis.line = element_line(size=2, colour = "black"),
      panel.grid.major = element_line(colour = "grey",size = 0.2,linetype = "dashed"),
      panel.grid.minor = element_line(colour = "grey",size = 0.2,linetype = "dashed"),
      panel.border = element_blank(), 
      panel.background = element_blank(),
      plot.title = element_text(size = 24,  face = "bold",hjust = 0.5),
      text=element_text(colour="black", size = 22),
      axis.text.x=element_text(colour="black", size = 22),
      axis.text.y=element_text(colour="black", size = 22),
      axis.title=element_text(size=24,colour = "black")) + 
    theme(strip.background =element_rect(fill="white")) + 
    theme(strip.text = element_text(size = 24))

# dev.off()

Table of means across conditions

# Make a table with conditions
posterior_means = as.data.frame(c("Attended Baseline High Reward", 
                                  "Attended Baseline Low Reward", 
                                  "Attended Acquisition High Reward", 
                                  "Attended Acquisition Low Reward", 
                                  "Attended Extinction High Reward", 
                                  "Attended Extinction Low Reward",
                                  "Unattended Baseline High Reward", 
                                  "Unattended Baseline Low Reward", 
                                  "Unattended Acquisition High Reward", 
                                  "Unattended Acquisition Low Reward", 
                                  "Unattended Extinction High Reward", 
                                  "Unattended Extinction Low Reward"))
names(posterior_means)[1] = "Condition"

posterior_means$Mean = c(paste(round(mean(Baseline_High_Attended), digits = 2), " [", round(hdi(Baseline_High_Attended)[[1]], digits = 2), " ", round(hdi(Baseline_High_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Baseline_Low_Attended), digits = 2), " [", round(hdi(Baseline_Low_Attended)[[1]], digits = 2), " ", round(hdi(Baseline_Low_Attended)[[2]], digits = 2), "]"),
                        
                        paste(round(mean(Acquisition_High_Attended), digits = 2), " [", round(hdi(Acquisition_High_Attended)[[1]], digits = 2), " ", round(hdi(Acquisition_High_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Acquisition_Low_Attended), digits = 2), " [", round(hdi(Acquisition_Low_Attended)[[1]], digits = 2), " ", round(hdi(Acquisition_Low_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_High_Attended), digits = 2), " [", round(hdi(Extinction_High_Attended)[[1]], digits = 2), " ", round(hdi(Extinction_High_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_Low_Attended), digits = 2), " [", round(hdi(Extinction_Low_Attended)[[1]], digits = 2), " ", round(hdi(Extinction_Low_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Baseline_High_NotAttended), digits = 2), " [", round(hdi(Baseline_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Baseline_High_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Baseline_Low_NotAttended), digits = 2), " [", round(hdi(Baseline_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Baseline_Low_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Acquisition_High_NotAttended), digits = 2), " [", round(hdi(Acquisition_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Acquisition_High_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Acquisition_Low_NotAttended), digits = 2), " [", round(hdi(Acquisition_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Acquisition_Low_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_High_NotAttended), digits = 2), " [", round(hdi(Extinction_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_High_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_Low_NotAttended), digits = 2), " [", round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2), "]"))

names(posterior_means)[2] = "Mean [HDI]"

posterior_means =  posterior_means %>% separate(Condition, c("Attention", "Reward Phase", "Reward Probability"), " ", extra = "merge")

kable(posterior_means, caption = "Means per condition")
Means per condition
Attention Reward Phase Reward Probability Mean [HDI]
Attended Baseline High Reward 1.11 [ 1.05 1.16 ]
Attended Baseline Low Reward 1.12 [ 1.07 1.17 ]
Attended Acquisition High Reward 1.14 [ 1.09 1.2 ]
Attended Acquisition Low Reward 1.1 [ 1.05 1.15 ]
Attended Extinction High Reward 1.11 [ 1.06 1.18 ]
Attended Extinction Low Reward 1.11 [ 1.05 1.18 ]
Unattended Baseline High Reward 0.86 [ 0.81 0.91 ]
Unattended Baseline Low Reward 0.87 [ 0.83 0.92 ]
Unattended Acquisition High Reward 0.9 [ 0.84 0.95 ]
Unattended Acquisition Low Reward 0.85 [ 0.81 0.9 ]
Unattended Extinction High Reward 0.87 [ 0.81 0.93 ]
Unattended Extinction Low Reward 0.87 [ 0.81 0.93 ]
### Inference about the best model

Attended vs. unattended

Check the difference between attended and not attended in baseline high rewarded

Diff_Att_NotAtt_Bsln_High = Baseline_High_Attended - Baseline_High_NotAttended
plotPost(Diff_Att_NotAtt_Bsln_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between attended and not attended in baseline low rewarded

Diff_Att_NotAtt_Bsln_Low = Baseline_Low_Attended - Baseline_Low_NotAttended
plotPost(Diff_Att_NotAtt_Bsln_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between attended and not attended in acquisition high rewarded

Diff_Att_NotAtt_Acq_High = Acquisition_High_Attended - Acquisition_High_NotAttended
plotPost(Diff_Att_NotAtt_Acq_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between attended and not attended in acquisition low rewarded

Diff_Att_NotAtt_Acq_Low = Acquisition_Low_Attended - Acquisition_Low_NotAttended
plotPost(Diff_Att_NotAtt_Acq_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between attended and not attended in extinction high rewarded

Diff_Att_NotAtt_Ext_High = Extinction_High_Attended - Extinction_High_NotAttended
plotPost(Diff_Att_NotAtt_Ext_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between attended and not attended in extinction low rewarded

Diff_Att_NotAtt_Ext_Low = Extinction_Low_Attended - Extinction_Low_NotAttended
plotPost(Diff_Att_NotAtt_Ext_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Comparison between phases

Check the difference between baseline and acquisition in high reward attended

Diff_Bsln_Acq_High_Att = Baseline_High_Attended - Acquisition_High_Attended
plotPost(Diff_Bsln_Acq_High_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between baseline and acquisition in low reward attended

Diff_Bsln_Acq_Low_Att = Baseline_Low_Attended - Acquisition_Low_Attended
plotPost(Diff_Bsln_Acq_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between baseline and acquisition in high reward not attended

Diff_Bsln_Acq_High_NotAtt = Baseline_High_NotAttended - Acquisition_High_NotAttended
plotPost(Diff_Bsln_Acq_High_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between baseline and acquisition in low reward not attended

Diff_Bsln_Acq_Low_NotAtt = Acquisition_Low_NotAttended - Baseline_Low_NotAttended  
plotPost(Diff_Bsln_Acq_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between acquisition and extinction in high reward attended

Diff_Acq_Ext_High_Att = Acquisition_High_Attended - Extinction_High_Attended
plotPost(Diff_Acq_Ext_High_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between acquisition and extinction in low reward attended

Diff_Acq_Ext_Low_Att = Extinction_Low_Attended - Acquisition_Low_Attended 
plotPost(Diff_Acq_Ext_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between acquisition and extinction in high reward not attended

Diff_Acq_Ext_High_NotAtt = Extinction_High_NotAttended - Acquisition_High_NotAttended 
plotPost(Diff_Acq_Ext_High_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between acquisition and extinction in low reward not attended

Diff_Acq_Ext_Low_NotAtt = Extinction_Low_NotAttended - Acquisition_Low_NotAttended 
plotPost(Diff_Acq_Ext_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Baseline difference

Check the difference between high and low reward in baseline attended

Diff_Bsln_High_Low_Att = Baseline_High_Attended - Baseline_Low_Attended
plotPost(Diff_Bsln_High_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between high and low reward in baseline not attended

Diff_Bsln_High_Low_NotAtt = Baseline_High_NotAttended - Baseline_Low_NotAttended
plotPost(Diff_Bsln_High_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between baseline and acquisition in high reward unattended vs. attended

Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended = Diff_Bsln_High_Low_NotAtt - Diff_Bsln_High_Low_Att 
plotPost(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

paste("Mean = ",round(mean(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended), digits = 2), " [", round(hdi(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended)[[1]], digits = 2), " ", round(hdi(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended)[[2]], digits = 2), "]")
## [1] "Mean =  0  [ 0   0 ]"

Only no movement trials

# Take only the movement trials
data = subset(data_all, Movement=="Nomov")

Means - raw data

summary = ddply(data,.(Attention,ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(Amplitude, na.rm = TRUE), digits = 2), " [", round(hdi(Amplitude)[[1]], digits = 2), " ", round(hdi(Amplitude)[[2]], digits = 2), "]")))

names(summary) = c("Attention", "Reward phase", "Reward probability", "Amplitude")

summary$Attention = recode(summary$Attention,
                           "Att" = "Attended",
                           "NotAtt" = "Unattended")

summary$`Reward phase` = recode(summary$`Reward phase`,
                           "Acq" = "Acquisition",
                           "Bsln" = "Baseline",
                           "Ext" = "Extinction")

summary$`Reward probability` = recode(summary$`Reward probability`,
                           "High_Rew" = "High",
                           "Low_Rew" = "Low")

summary = as.data.frame(summary)

summary$`Reward phase` = factor(summary$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
summary = summary[order(summary$Attention,summary$`Reward phase`,summary$`Reward probability`),]
row.names(summary) = NULL

kable(summary, caption = "Amplitudes per condition")
Amplitudes per condition
Attention Reward phase Reward probability Amplitude
Attended Baseline High 1.16 [ 0.99 1.52 ]
Attended Baseline Low 1.14 [ 0.9 1.44 ]
Attended Acquisition High 1.14 [ 0.8 1.49 ]
Attended Acquisition Low 1.11 [ 0.76 1.56 ]
Attended Extinction High 1.11 [ 0.64 1.43 ]
Attended Extinction Low 1.14 [ 0.85 1.97 ]
Unattended Baseline High 0.87 [ 0.62 1.22 ]
Unattended Baseline Low 0.87 [ 0.6 1.27 ]
Unattended Acquisition High 0.91 [ 0.54 1.36 ]
Unattended Acquisition Low 0.89 [ 0.62 1.43 ]
Unattended Extinction High 0.88 [ 0.55 1.23 ]
Unattended Extinction Low 0.91 [ 0.5 1.42 ]

Plots - raw data

# Plot amplitude across experiment phases------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# prepare data for plotting
dataPlot = data

# rename variables
colnames(dataPlot)[colnames(dataPlot)=="ExpPhase"] <- "Reward phase"
colnames(dataPlot)[colnames(dataPlot)=="Condition"] <- "Reward probability"

# rename conditions
dataPlot$`Reward phase` = recode(dataPlot$`Reward phase`,
                                  "Acq" = "Acquisition",
                                  "Bsln" = "Baseline",
                                  "Ext" = "Extinction")

dataPlot$`Reward probability` = recode(dataPlot$`Reward probability`,
                                        "High_Rew" = "High",
                                        "Low_Rew" = "Low")

#order
dataPlot$`Reward phase` = factor(dataPlot$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward phase`,dataPlot$`Reward probability`),]


plottingConditions = c("Attended","Unattended" )
for (i in 1:length(plottingConditions)){
  
  if(plottingConditions[i]=="Attended"){dataAmplitudePlot=subset(dataPlot,Attention=="Att")}
  
  if(plottingConditions[i]=="Unattended"){dataAmplitudePlot=subset(dataPlot,Attention=="NotAtt")}  

# Pirate plot

    pirateplot(formula = Amplitude ~ `Reward phase` + `Reward probability`, # dependent~independent variables
             data=dataAmplitudePlot, # data frame
             main=plottingConditions[i], # main title
             ylim=c(0.2,2.2), # y-axis: limits
             ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
             theme=0, # preset theme (0: use your own)
             point.col="black", # points: color
             point.o=.3, # points: opacity (0-1)
             avg.line.col="black", # average line: color
             avg.line.lwd=2, # average line: line width
             avg.line.o=1, # average line: opacity (0-1)
             bean.b.col="black", # bean border, color
             bean.lwd=0.6, # bean border, line width
             bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
             bean.b.o=0.3, # bean border, opacity (0-1)
             bean.f.col="gray", # bean filling, color
             bean.f.o=.1, # bean filling, opacity (0-1)
             cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
             gl.col="gray", # gridlines: color
             gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
             cex.lab=1, # axis labels: size
             cex.axis=1, # axis numbers: size
             cex.names = 1,
             bty="l", # plot box type
             back.col="white") # background, color
}

Statistics

# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))

model.rewardTimesPhasePlusAtt = readRDS("rewardTimesPhasePlusAtt.EEG.nomovementtrials.rds")
compare.threefactors.EEG.waic = readRDS("compare.EEG.waic.nomovementtrials.rds")
bR2.null.EEG = readRDS("bR2.null.EEG.nomovementtrials")
bR2.expphase.EEG = readRDS("bR2.expphase.EEG.nomovementtrials")
bR2.attention.EEG = readRDS("bR2.attention.EEG.nomovementtrials")
bR2.phaseANDattention.EEG = readRDS("bR2.phaseANDattention.EEG.nomovementtrials")
bR2.phaseANDattention_interaction.EEG = readRDS("bR2.phaseANDattention_interaction.EEG.nomovementtrials")
bR2.rewardTimesPhasePlusAtt.EEG = readRDS("bR2.rewardTimesPhasePlusAtt.EEG.nomovementtrials")
bR2.full.EEG = readRDS("bR2.full.EEG.nomovementtrials")

Model comparison with WAIC

print(compare.threefactors.EEG.waic)
## Output of model 'null':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic     -7.0 20.0
## p_waic        17.1  1.5
## waic          13.9 39.9
## 
## Output of model 'expphase':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic    -10.8 20.0
## p_waic        25.2  2.4
## waic          21.6 40.0
## 
## Output of model 'attention':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic     91.6 22.8
## p_waic        53.4  4.6
## waic        -183.2 45.5
## 
## Output of model 'phaseANDattention':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic     87.8 22.8
## p_waic        65.0  6.1
## waic        -175.6 45.7
## 
## Output of model 'phaseANDattention_interaction':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic     86.7 22.9
## p_waic        66.8  6.2
## waic        -173.4 45.8
## 
## Output of model 'rewardTimesPhasePlusAtt':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic    128.0 24.7
## p_waic       102.7  9.5
## waic        -256.1 49.5
## 
## Output of model 'full':
## 
## Computed from 16000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic    123.7 24.8
## p_waic       106.7  9.9
## waic        -247.4 49.6
## 
## Model comparisons:
##                               elpd_diff se_diff
## rewardTimesPhasePlusAtt          0.0       0.0 
## full                            -4.3       1.6 
## attention                      -36.4      11.2 
## phaseANDattention              -40.2      11.2 
## phaseANDattention_interaction  -41.3      11.4 
## null                          -135.0      16.4 
## expphase                      -138.8      16.5

Bayesian R squared

Null model

print(bR2.null.EEG)
##      Estimate  Est.Error         Q2.5      Q97.5
## R2 0.03957342 0.02418746 0.0007912515 0.09155692

Experiment phase

print(bR2.expphase.EEG)
##      Estimate Est.Error       Q2.5     Q97.5
## R2 0.05514417 0.0252663 0.01162466 0.1079302

Attention model

print(bR2.attention.EEG)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.3952447 0.03200049 0.3281384 0.4539408

Phase and attention model

print(bR2.phaseANDattention.EEG)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.4079176 0.03142196 0.3419071 0.4650304

Phase and attention interaction model

print(bR2.phaseANDattention_interaction.EEG)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.4106584 0.03145109 0.3451474 0.4687237

Reward times phase plus attention model

print(bR2.rewardTimesPhasePlusAtt.EEG)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.5403708 0.02781604 0.4815786 0.5899214

Full model

print(bR2.full.EEG)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.5414225 0.02735935 0.4838342 0.5907724

Checking the best model

Plotting the chains

# Plot chains
plot(model.rewardTimesPhasePlusAtt, pars = "^b_", ask = FALSE, N=6)

Summary of the best model

# Summary of the best model
print(tidy(model.rewardTimesPhasePlusAtt, par_type = "non-varying"), digits = 1)
##                           term estimate std.error lower upper
## 1                    Intercept    1.137      0.02  1.10  1.18
## 2             ConditionLow_Rew   -0.009      0.04 -0.07  0.05
## 3                  ExpPhaseAcq    0.012      0.03 -0.03  0.05
## 4                  ExpPhaseExt   -0.023      0.03 -0.07  0.02
## 5              AttentionNotAtt   -0.243      0.03 -0.29 -0.20
## 6 ConditionLow_Rew:ExpPhaseAcq   -0.014      0.04 -0.07  0.04
## 7 ConditionLow_Rew:ExpPhaseExt    0.040      0.04 -0.02  0.10

Plotting the best model

post = posterior_samples(model.rewardTimesPhasePlusAtt, "^b")

# Calculate posteriors for each condition

################################################ Baseline ####

##################### Attended

######### High reward
Baseline_High_Attended = post[["b_Intercept"]]
######### Low reward
Baseline_Low_Attended = post[["b_Intercept"]] + 
  post[["b_ConditionLow_Rew"]] 

##################### Not Attended

######### High reward
Baseline_High_NotAttended = post[["b_Intercept"]] + 
  post[["b_AttentionNotAtt"]]
######### Low reward
Baseline_Low_NotAttended = post[["b_Intercept"]] + 
  post[["b_AttentionNotAtt"]] + 
  post[["b_ConditionLow_Rew"]] 

################################################ Acquistion

##################### Attended

######### High reward
Acquisition_High_Attended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] 
######### Low reward
Acquisition_Low_Attended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ConditionLow_Rew:ExpPhaseAcq"]]

##################### Not Attended

######### High reward
Acquisition_High_NotAttended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] + 
  post[["b_AttentionNotAtt"]] 

######### Low reward
Acquisition_Low_NotAttended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] + 
  post[["b_AttentionNotAtt"]] + 
  post[["b_ConditionLow_Rew"]] +
  post[["b_ConditionLow_Rew:ExpPhaseAcq"]] 


################################################ Extinction

##################### Attended

######### High reward
Extinction_High_Attended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] 
######### Low reward
Extinction_Low_Attended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ConditionLow_Rew:ExpPhaseExt"]]

##################### Not Attended

######### High reward
Extinction_High_NotAttended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] + 
  post[["b_AttentionNotAtt"]]
######### Low reward
Extinction_Low_NotAttended = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] + 
  post[["b_AttentionNotAtt"]] + 
  post[["b_ConditionLow_Rew"]] +
  post[["b_ConditionLow_Rew:ExpPhaseExt"]] 
# make a data frame

posterior_conditions = melt(data.frame(Baseline_High_Attended, Baseline_High_NotAttended, Baseline_Low_Attended, Baseline_Low_NotAttended, Acquisition_High_Attended, Acquisition_High_NotAttended, Acquisition_Low_Attended, Acquisition_Low_NotAttended, Extinction_High_Attended, Extinction_High_NotAttended, Extinction_Low_Attended, Extinction_Low_NotAttended))

posterior_conditions =  posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability", "Attention"), "_", extra = "merge")

posterior_conditions$Attention = recode(posterior_conditions$Attention,
                           "Attended" = "Attended",
                           "NotAttended" = "Unattended")

names(posterior_conditions)[4] = "Amplitude"


#order
#dataPlot$`Reward phase` = factor(dataPlot$`Reward phase`, levels = c("Baseline","Acquisition","Extinction"))
#dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward phase`,dataPlot$`Reward probability`),]


plottingConditions = c("Attended","Unattended" )
for (i in 1:length(plottingConditions)){
  
  if(plottingConditions[i]=="Attended"){dataAmplitudePlot=subset(posterior_conditions,Attention=="Attended")}
  
  if(plottingConditions[i]=="Unattended"){dataAmplitudePlot=subset(posterior_conditions,Attention=="Unattended")}  

# Pirate plot

    pirateplot(formula = Amplitude ~ `Reward Phase` + `Reward Probability`, # dependent~independent variables
             data=dataAmplitudePlot, # data frame
             main=plottingConditions[i], # main title
             ylim=c(0.7,1.3), # y-axis: limits
             ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
             theme=0, # preset theme (0: use your own)
             avg.line.col="black", # average line: color
             avg.line.lwd=2, # average line: line width
             avg.line.o=1, # average line: opacity (0-1)
             bean.b.col="black", # bean border, color
             bean.lwd=0.6, # bean border, line width
             bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
             bean.b.o=0.3, # bean border, opacity (0-1)
             bean.f.col="gray", # bean filling, color
             bean.f.o=.1, # bean filling, opacity (0-1)
             cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
             gl.col="gray", # gridlines: color
             gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
             cex.lab=1, # axis labels: size
             cex.axis=1, # axis numbers: size
             cex.names = 1,
             sortx = "sequential",
             bty="l", # plot box type
             back.col="white") # background, color
}

plottingConditions = c("High","Low" )
for (i in 1:length(plottingConditions)){
  
  if(plottingConditions[i]=="High"){dataAmplitudePlot=subset(posterior_conditions,`Reward Probability`=="High")}
  
  if(plottingConditions[i]=="Low"){dataAmplitudePlot=subset(posterior_conditions,`Reward Probability`=="Low")}  

# Pirate plot

    pirateplot(formula = Amplitude ~ `Attention` + `Reward Phase`, # dependent~independent variables
             data=dataAmplitudePlot, # data frame
             main=plottingConditions[i], # main title
             ylim=c(0.7,1.3), # y-axis: limits
             ylab=expression(paste("Amplitude (",mu,"V)")), # y-axis: label
             theme=0, # preset theme (0: use your own)
             avg.line.col="black", # average line: color
             avg.line.lwd=2, # average line: line width
             avg.line.o=1, # average line: opacity (0-1)
             bean.b.col="black", # bean border, color
             bean.lwd=0.6, # bean border, line width
             bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
             bean.b.o=0.3, # bean border, opacity (0-1)
             bean.f.col="gray", # bean filling, color
             bean.f.o=.1, # bean filling, opacity (0-1)
             cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
             gl.col="gray", # gridlines: color
             gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
             cex.lab=1, # axis labels: size
             cex.axis=1, # axis numbers: size
             cex.names = 1,
             sortx = "sequential",
             bty="l", # plot box type
             back.col="white") # background, color
}

# Data from the model

# make a data frame of conditions from the posterior

posterior_conditions = melt(data.frame(Baseline_High_Attended, Baseline_High_NotAttended, Baseline_Low_Attended, Baseline_Low_NotAttended, Acquisition_High_Attended, Acquisition_High_NotAttended, Acquisition_Low_Attended, Acquisition_Low_NotAttended, Extinction_High_Attended, Extinction_High_NotAttended, Extinction_Low_Attended, Extinction_Low_NotAttended))

posterior_conditions =  posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability", "Attention"), "_", extra = "merge")

posterior_conditions$Attention = recode(posterior_conditions$Attention,
                           "Attended" = "Target",
                           "NotAttended" = "Distractor")

posterior_conditions$Attention  = as.factor(posterior_conditions$Attention )

posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
                                      "High" = "High reward probability",
                                      "Low" = "Low reward probability")

posterior_conditions$`Reward Probability` = as.factor(posterior_conditions$`Reward Probability`)

posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
                                "Acquisition" = "Training",
                                "Baseline" = "Baseline",
                                "Extinction" = "Test")

posterior_conditions$`Reward Phase` = as.factor(posterior_conditions$`Reward Phase`)

names(posterior_conditions)[4] = "Amplitude"




# Make a table with credible intervals
posterior_means = as.data.frame(c("Attended Baseline High Reward", 
                                  "Unattended Baseline High Reward", 
                                  "Attended Baseline Low Reward",
                                  "Unattended Baseline Low Reward",
                                  
                                  "Attended Acquisition High Reward", 
                                  "Unattended Acquisition High Reward", 
                                  "Attended Acquisition Low Reward",
                                  "Unattended Acquisition Low Reward",
                                  
                                  "Attended Extinction High Reward", 
                                  "Unattended Extinction High Reward", 
                                  "Attended Extinction Low Reward",
                                  "Unattended Extinction Low Reward"))

names(posterior_means)[1] = "Condition"

posterior_means$low_95_CI = c(round(hdi(Baseline_High_Attended)[[1]], digits = 2),
                              round(hdi(Baseline_High_NotAttended)[[1]], digits = 2),
                              round(hdi(Baseline_Low_Attended)[[1]], digits = 2),
                              round(hdi(Baseline_Low_NotAttended)[[1]], digits = 2),
                              round(hdi(Acquisition_High_Attended)[[1]], digits = 2),
                              round(hdi(Acquisition_High_NotAttended)[[1]], digits = 2),
                              round(hdi(Acquisition_Low_Attended)[[1]], digits = 2),
                              round(hdi(Acquisition_Low_NotAttended)[[1]], digits = 2),
                              round(hdi(Extinction_High_Attended)[[1]], digits = 2),
                              round(hdi(Extinction_High_NotAttended)[[1]], digits = 2),
                              round(hdi(Extinction_Low_Attended)[[1]], digits = 2),
                              round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2))

posterior_means$high_95_CI = c(round(hdi(Baseline_High_Attended)[[2]], digits = 2),
                              round(hdi(Baseline_High_NotAttended)[[2]], digits = 2),
                              round(hdi(Baseline_Low_Attended)[[2]], digits = 2),
                              round(hdi(Baseline_Low_NotAttended)[[2]], digits = 2),
                              round(hdi(Acquisition_High_Attended)[[2]], digits = 2),
                              round(hdi(Acquisition_High_NotAttended)[[2]], digits = 2),
                              round(hdi(Acquisition_Low_Attended)[[2]], digits = 2),
                              round(hdi(Acquisition_Low_NotAttended)[[2]], digits = 2),
                              round(hdi(Extinction_High_Attended)[[2]], digits = 2),
                              round(hdi(Extinction_High_NotAttended)[[2]], digits = 2),
                              round(hdi(Extinction_Low_Attended)[[2]], digits = 2),
                              round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2))

posterior_means$low_50_CI = c(round(hdi(Baseline_High_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Baseline_High_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Baseline_Low_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Baseline_Low_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Acquisition_High_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Acquisition_High_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Acquisition_Low_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Acquisition_Low_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Extinction_High_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Extinction_High_NotAttended,0.5)[[1]], digits = 2),
                              round(hdi(Extinction_Low_Attended,0.5)[[1]], digits = 2),
                              round(hdi(Extinction_Low_NotAttended,0.5)[[1]], digits = 2))

posterior_means$high_50_CI = c(round(hdi(Baseline_High_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Baseline_High_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Baseline_Low_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Baseline_Low_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Acquisition_High_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Acquisition_High_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Acquisition_Low_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Acquisition_Low_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Extinction_High_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Extinction_High_NotAttended,0.5)[[2]], digits = 2),
                              round(hdi(Extinction_Low_Attended,0.5)[[2]], digits = 2),
                              round(hdi(Extinction_Low_NotAttended,0.5)[[2]], digits = 2))


posterior_means =  posterior_means %>% separate(Condition, c("Attention","Reward Phase", "Reward Probability"), " ", extra = "merge")

posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
                                      "High Reward" = "High reward probability",
                                      "Low Reward" = "Low reward probability")

posterior_means$`Reward Probability`= as.factor(posterior_means$`Reward Probability`)

posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
                                "Acquisition" = "Training",
                                "Baseline" = "Baseline",
                                "Extinction" = "Test")

posterior_means$`Reward Phase` = as.factor(posterior_means$`Reward Phase`)

posterior_means$Attention = recode(posterior_means$Attention,
                           "Attended" = "Target",
                           "Unattended" = "Distractor")

posterior_means$Attention = as.factor(posterior_means$Attention)

# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$Amplitude = c(1,1,1,1,1,1,1,1,1,1,1,1)


# Raw data
# Prepare the dataset
# prepare data for plotting - average over motion

dataPlot = data_all
dataPlot = ddply(dataPlot,.(Attention,ExpPhase,Condition,Subject),plyr::summarize,Amplitude=mean(Amplitude, na.rm = TRUE))


# rename variables
colnames(dataPlot)[colnames(dataPlot)=="ExpPhase"] <- "Reward Phase"
colnames(dataPlot)[colnames(dataPlot)=="Condition"] <- "Reward Probability"

# rename conditions
dataPlot$`Reward Phase` = recode(dataPlot$`Reward Phase`,
                                  "Acq" = "Training",
                                  "Bsln" = "Baseline",
                                  "Ext" = "Test")

dataPlot$`Reward Probability` = recode(dataPlot$`Reward Probability`,
                                        "High_Rew" = "High reward probability",
                                        "Low_Rew" = "Low reward probability")

dataPlot$Attention = recode(dataPlot$Attention,
                                        "Att" = "Target",
                                        "NotAtt" = "Distractor")


##### Plotting targets and distractors separately  #####
  
# Attended
  
  dataAmplitudePlot=subset(dataPlot,Attention=="Target")
  posterior_means_plot=subset(posterior_means,Attention=="Target")
  posterior_conditions_plot=subset(posterior_conditions,Attention=="Target")
  
  # tiff("ssvep_attended.tiff", units="in", width=8, height=5, res=800)
    
  ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) + 
    # violin of the raw data
    geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
    # individual participants for the raw data
    # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
    geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
    # mean of the model
    stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
    # 95 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
    # 50 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
    scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
    facet_grid(.~`Reward Probability`) + 
    scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
    theme(
      axis.line = element_line(size=1, colour = "black"),
      panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.border = element_blank(), 
      panel.background = element_blank(),
      plot.title = element_text(size = 12,  face = "bold",hjust = 0.5),
      text=element_text(colour="black", size = 11),
      axis.text.x=element_text(colour="black", size = 11),
      axis.text.y=element_text(colour="black", size = 11),
      axis.title=element_text(size=12,colour = "black")) + 
    theme(strip.background =element_rect(fill="white")) + 
    theme(strip.text = element_text(size = 12))

  # dev.off()
  
  
  
  # Unattended
  
  dataAmplitudePlot=subset(dataPlot,Attention=="Distractor")
  posterior_means_plot=subset(posterior_means,Attention=="Distractor")
  posterior_conditions_plot=subset(posterior_conditions,Attention=="Distractor")
  
  # tiff("ssvep_unattended.tiff", units="in", width=8, height=5, res=800)
    
  ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) + 
    # violin of the raw data
    geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
    # individual participants for the raw data
    # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
    geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
    # mean of the model
    stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
    # 95 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
    # 50 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
    scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
    facet_grid(.~`Reward Probability`) + 
    scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
    theme(
      axis.line = element_line(size=1, colour = "black"),
      panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.border = element_blank(), 
      panel.background = element_blank(),
      plot.title = element_text(size = 12,  face = "bold",hjust = 0.5),
      text=element_text(colour="black", size = 11),
      axis.text.x=element_text(colour="black", size = 11),
      axis.text.y=element_text(colour="black", size = 11),
      axis.title=element_text(size=12,colour = "black")) + 
    theme(strip.background =element_rect(fill="white")) + 
    theme(strip.text = element_text(size = 12))

  # dev.off()

  
  
##### Plotting high and low separately  #####
  
# High
  
  dataAmplitudePlot=subset(dataPlot,`Reward Probability`=="High reward probability")
  posterior_means_plot=subset(posterior_means,`Reward Probability`=="High reward probability")
  posterior_conditions_plot=subset(posterior_conditions,`Reward Probability`=="High reward probability")
  
  dataPlot$Attention = relevel(dataPlot$Attention,"Target")
  posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
  posterior_means$Attention = relevel(posterior_means$Attention,"Target")
  
  # tiff("ssvep_high.tiff", units="in", width=8, height=5, res=800)
    
  ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=Attention)) + 
    # violin of the raw data
    geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
    # individual participants for the raw data
    # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
    geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
    # mean of the model
    stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=Attention),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
    # 95 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
    # 50 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
    scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
    facet_grid(.~Attention) + 
    scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
    labs(title="High reward") +
    theme(
      axis.line = element_line(size=1, colour = "black"),
      panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.border = element_blank(), 
      panel.background = element_blank(),
      plot.title = element_text(size = 12,  face = "bold",hjust = 0.5),
      text=element_text(colour="black", size = 11),
      axis.text.x=element_text(colour="black", size = 11),
      axis.text.y=element_text(colour="black", size = 11),
      axis.title=element_text(size=12,colour = "black")) + 
    theme(strip.background =element_rect(fill="white")) + 
    theme(strip.text = element_text(size = 12))

  # dev.off()
  
  
  
# Low
  
  dataAmplitudePlot=subset(dataPlot,`Reward Probability`=="Low reward probability")
  posterior_means_plot=subset(posterior_means,`Reward Probability`=="Low reward probability")
  posterior_conditions_plot=subset(posterior_conditions,`Reward Probability`=="Low reward probability")
  
  dataPlot$Attention = relevel(dataPlot$Attention,"Target")
  posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
  posterior_means$Attention = relevel(posterior_means$Attention,"Target")
  
  # tiff("ssvep_low.tiff", units="in", width=8, height=5, res=800)
    
  ggplot(data=dataAmplitudePlot, aes(x=`Reward Phase`, y=`Amplitude`, label=Attention)) + 
    # violin of the raw data
    geom_violin(data=dataAmplitudePlot, color ="#636363",trim = TRUE,alpha="0.2") +
    # individual participants for the raw data
    # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
    geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
    # mean of the model
    stat_summary(data = posterior_conditions_plot, aes(x=`Reward Phase`, y=`Amplitude`, group=Attention),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
    # 95 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
    # 50 credible intervals of the model
    geom_linerange(data = posterior_means_plot,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
    scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
    facet_grid(.~Attention) + 
    scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(-1,3,0.5)) +
    labs(title="Low reward") +
    theme(
      axis.line = element_line(size=1, colour = "black"),
      panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
      panel.border = element_blank(), 
      panel.background = element_blank(),
      plot.title = element_text(size = 12,  face = "bold",hjust = 0.5),
      text=element_text(colour="black", size = 11),
      axis.text.x=element_text(colour="black", size = 11),
      axis.text.y=element_text(colour="black", size = 11),
      axis.title=element_text(size=12,colour = "black")) + 
    theme(strip.background =element_rect(fill="white")) + 
    theme(strip.text = element_text(size = 12))

  # dev.off()
####### Plotting everything together #####
  
# Order for plotting
# dataPlot$Attention = factor(dataPlot$Attention, levels = c("Target","Distractor")) 
# posterior_conditions$Attention = factor(posterior_conditions$Attention, levels = c("Target","Distractor")) 
# posterior_means$Attention = factor(posterior_means$Attention, levels = c("Target","Distractor")) 
# 
# dataPlot = dataPlot[order(dataPlot$Attention,dataPlot$`Reward Phase`,dataPlot$`Reward Probability`),]

dataPlot$Attention = relevel(dataPlot$Attention,"Target")
posterior_conditions$Attention = relevel(posterior_conditions$Attention,"Target")
posterior_means$Attention = relevel(posterior_means$Attention,"Target")


  
  
    # tiff("ssvep_both.tiff", units="in", width=16, height=10, res=800)

    ggplot(data=dataPlot, aes(x=`Reward Phase`, y=`Amplitude`, label=`Reward Probability`)) + 
    # violin of the raw data
    geom_violin(data=dataPlot, color ="#636363",trim = TRUE,alpha="0.2") +
    # individual participants for the raw data
    # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
    geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
    # mean of the model
    stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`Amplitude`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
    # 95 credible intervals of the model
    geom_linerange(data = posterior_means,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_95_CI, ymax=`Amplitude`-`Amplitude`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
    # 50 credible intervals of the model
    geom_linerange(data = posterior_means,width=.1, aes(ymin=`Amplitude`-`Amplitude`+low_50_CI, ymax=`Amplitude`-`Amplitude`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
    scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
    facet_grid(`Reward Probability`~Attention) + 
    scale_y_continuous(name="Amplitude (a.u.)",breaks=seq(0.5,2,0.5)) +
    theme(
      axis.line = element_line(size=2, colour = "black"),
      panel.grid.major = element_line(colour = "grey",size = 0.2,linetype = "dashed"),
      panel.grid.minor = element_line(colour = "grey",size = 0.2,linetype = "dashed"),
      panel.border = element_blank(), 
      panel.background = element_blank(),
      plot.title = element_text(size = 24,  face = "bold",hjust = 0.5),
      text=element_text(colour="black", size = 22),
      axis.text.x=element_text(colour="black", size = 22),
      axis.text.y=element_text(colour="black", size = 22),
      axis.title=element_text(size=24,colour = "black")) + 
    theme(strip.background =element_rect(fill="white")) + 
    theme(strip.text = element_text(size = 24))

# dev.off()

Table of means across conditions

# Make a table with conditions
posterior_means = as.data.frame(c("Attended Baseline High Reward", 
                                  "Attended Baseline Low Reward", 
                                  "Attended Acquisition High Reward", 
                                  "Attended Acquisition Low Reward", 
                                  "Attended Extinction High Reward", 
                                  "Attended Extinction Low Reward",
                                  "Unattended Baseline High Reward", 
                                  "Unattended Baseline Low Reward", 
                                  "Unattended Acquisition High Reward", 
                                  "Unattended Acquisition Low Reward", 
                                  "Unattended Extinction High Reward", 
                                  "Unattended Extinction Low Reward"))
names(posterior_means)[1] = "Condition"

posterior_means$Mean = c(paste(round(mean(Baseline_High_Attended), digits = 2), " [", round(hdi(Baseline_High_Attended)[[1]], digits = 2), " ", round(hdi(Baseline_High_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Baseline_Low_Attended), digits = 2), " [", round(hdi(Baseline_Low_Attended)[[1]], digits = 2), " ", round(hdi(Baseline_Low_Attended)[[2]], digits = 2), "]"),
                        
                        paste(round(mean(Acquisition_High_Attended), digits = 2), " [", round(hdi(Acquisition_High_Attended)[[1]], digits = 2), " ", round(hdi(Acquisition_High_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Acquisition_Low_Attended), digits = 2), " [", round(hdi(Acquisition_Low_Attended)[[1]], digits = 2), " ", round(hdi(Acquisition_Low_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_High_Attended), digits = 2), " [", round(hdi(Extinction_High_Attended)[[1]], digits = 2), " ", round(hdi(Extinction_High_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_Low_Attended), digits = 2), " [", round(hdi(Extinction_Low_Attended)[[1]], digits = 2), " ", round(hdi(Extinction_Low_Attended)[[2]], digits = 2), "]"),
                        paste(round(mean(Baseline_High_NotAttended), digits = 2), " [", round(hdi(Baseline_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Baseline_High_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Baseline_Low_NotAttended), digits = 2), " [", round(hdi(Baseline_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Baseline_Low_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Acquisition_High_NotAttended), digits = 2), " [", round(hdi(Acquisition_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Acquisition_High_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Acquisition_Low_NotAttended), digits = 2), " [", round(hdi(Acquisition_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Acquisition_Low_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_High_NotAttended), digits = 2), " [", round(hdi(Extinction_High_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_High_NotAttended)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_Low_NotAttended), digits = 2), " [", round(hdi(Extinction_Low_NotAttended)[[1]], digits = 2), " ", round(hdi(Extinction_Low_NotAttended)[[2]], digits = 2), "]"))

names(posterior_means)[2] = "Mean [HDI]"

posterior_means =  posterior_means %>% separate(Condition, c("Attention", "Reward Phase", "Reward Probability"), " ", extra = "merge")

kable(posterior_means, caption = "Means per condition")
Means per condition
Attention Reward Phase Reward Probability Mean [HDI]
Attended Baseline High Reward 1.14 [ 1.09 1.18 ]
Attended Baseline Low Reward 1.13 [ 1.07 1.18 ]
Attended Acquisition High Reward 1.15 [ 1.1 1.2 ]
Attended Acquisition Low Reward 1.13 [ 1.07 1.19 ]
Attended Extinction High Reward 1.11 [ 1.06 1.17 ]
Attended Extinction Low Reward 1.14 [ 1.08 1.21 ]
Unattended Baseline High Reward 0.89 [ 0.84 0.95 ]
Unattended Baseline Low Reward 0.88 [ 0.82 0.94 ]
Unattended Acquisition High Reward 0.91 [ 0.85 0.96 ]
Unattended Acquisition Low Reward 0.88 [ 0.83 0.94 ]
Unattended Extinction High Reward 0.87 [ 0.81 0.93 ]
Unattended Extinction Low Reward 0.9 [ 0.84 0.97 ]
### Inference about the best model

Attended vs. unattended

Check the difference between attended and not attended in baseline high rewarded

Diff_Att_NotAtt_Bsln_High = Baseline_High_Attended - Baseline_High_NotAttended
plotPost(Diff_Att_NotAtt_Bsln_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between attended and not attended in baseline low rewarded

Diff_Att_NotAtt_Bsln_Low = Baseline_Low_Attended - Baseline_Low_NotAttended
plotPost(Diff_Att_NotAtt_Bsln_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between attended and not attended in acquisition high rewarded

Diff_Att_NotAtt_Acq_High = Acquisition_High_Attended - Acquisition_High_NotAttended
plotPost(Diff_Att_NotAtt_Acq_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between attended and not attended in acquisition low rewarded

Diff_Att_NotAtt_Acq_Low = Acquisition_Low_Attended - Acquisition_Low_NotAttended
plotPost(Diff_Att_NotAtt_Acq_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between attended and not attended in extinction high rewarded

Diff_Att_NotAtt_Ext_High = Extinction_High_Attended - Extinction_High_NotAttended
plotPost(Diff_Att_NotAtt_Ext_High, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between attended and not attended in extinction low rewarded

Diff_Att_NotAtt_Ext_Low = Extinction_Low_Attended - Extinction_Low_NotAttended
plotPost(Diff_Att_NotAtt_Ext_Low, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Comparison between phases

Check the difference between baseline and acquisition in high reward attended

Diff_Bsln_Acq_High_Att = Baseline_High_Attended - Acquisition_High_Attended
plotPost(Diff_Bsln_Acq_High_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between baseline and acquisition in low reward attended

Diff_Bsln_Acq_Low_Att = Baseline_Low_Attended - Acquisition_Low_Attended
plotPost(Diff_Bsln_Acq_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between baseline and acquisition in high reward not attended

Diff_Bsln_Acq_High_NotAtt = Baseline_High_NotAttended - Acquisition_High_NotAttended
plotPost(Diff_Bsln_Acq_High_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between baseline and acquisition in low reward not attended

Diff_Bsln_Acq_Low_NotAtt = Acquisition_Low_NotAttended - Baseline_Low_NotAttended  
plotPost(Diff_Bsln_Acq_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between acquisition and extinction in high reward attended

Diff_Acq_Ext_High_Att = Acquisition_High_Attended - Extinction_High_Attended
plotPost(Diff_Acq_Ext_High_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between acquisition and extinction in low reward attended

Diff_Acq_Ext_Low_Att = Extinction_Low_Attended - Acquisition_Low_Attended 
plotPost(Diff_Acq_Ext_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between acquisition and extinction in high reward not attended

Diff_Acq_Ext_High_NotAtt = Extinction_High_NotAttended - Acquisition_High_NotAttended 
plotPost(Diff_Acq_Ext_High_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between acquisition and extinction in low reward not attended

Diff_Acq_Ext_Low_NotAtt = Extinction_Low_NotAttended - Acquisition_Low_NotAttended 
plotPost(Diff_Acq_Ext_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Baseline difference

Check the difference between high and low reward in baseline attended

Diff_Bsln_High_Low_Att = Baseline_High_Attended - Baseline_Low_Attended
plotPost(Diff_Bsln_High_Low_Att, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between high and low reward in baseline not attended

Diff_Bsln_High_Low_NotAtt = Baseline_High_NotAttended - Baseline_Low_NotAttended
plotPost(Diff_Bsln_High_Low_NotAtt, xlab = "", col = "#b3cde0", cex = 1, showCurve = FALSE, compVal = 0)

Check the difference between baseline and acquisition in high reward unattended vs. attended

Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended = Diff_Bsln_High_Low_NotAtt - Diff_Bsln_High_Low_Att 
plotPost(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

paste("Mean = ",round(mean(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended), digits = 2), " [", round(hdi(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended)[[1]], digits = 2), " ", round(hdi(Diff_Bsln_Acq_High_vs_Low_in_attended_vs_unattended)[[2]], digits = 2), "]")
## [1] "Mean =  0  [ 0   0 ]"

Behavior

Import data

# Clear environemnt and import data------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

# # clear the environment
# rm(list=ls()) 
# # clear the console
# cat("\014") 
# #load packages and install them if they're not installed
# if (!require("pacman")) install.packages("pacman")
# pacman::p_load(plyr,Rmisc,yarrr,BayesFactor,reshape2,brms, broom, tidyverse, brmstools, BEST, knitr, here,psych)
# # set seed
# set.seed(42) 
# # set directory
# setwd(here())
# import data
data.raw = read.csv(file = here("data","Data_behavior_exp1_48pps.csv"),header=TRUE,na.strings="NaN")

# Prepare the dataset------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

### Adding and renaming variables 
# rename EventType variable
names(data.raw)[names(data.raw) == "EventType"] = "MovedDots" 
# add a variable with the name of the attended color instead of a numbers
data.raw$AttendedColor = ifelse(data.raw$AttendedColor==1,"red","blue")
# add a variable saying which color was linked with High_Rew (even numbers - blue was High_Rew)
data.raw$RewardedColor = ifelse(data.raw$ParticipantNo%%2==0,"blue","red") 
# add a variable with the name of the moved color instead of a numbers
data.raw$MovedDots = ifelse(data.raw$MovedDots==1,"red","blue") 
# split experimental phases into 6 isntead of 3 phases (trial 0-200: Bsln; trial 201-400: Acq; trial 401-600: Ext)
#data.raw$ExpPhase = cut(data.raw$Trial,breaks=c(0,100,200,300,400,500,600),labels=c("Bsln1","Bsln2","Acq1","Acq2","Ext1","Ext2")) 
# split experimental phases into 3 phases (trial 0-200: Bsln; trial 201-400: Acq; trial 401-600: Ext)
data.raw$ExpPhase = cut(data.raw$Trial,breaks=c(0,200,400,600),labels=c("Bsln","Acq","Ext")) # trial 0-200: Bsln; trial 201-400: Acq; trial 401-600: Ext

### Convert variables to be used in analyses into factors
data.raw[c("ParticipantNo", "AttendedColor","RewardedColor", "MovedDots", "ExpPhase" )] = 
  lapply(data.raw[c("ParticipantNo", "AttendedColor","RewardedColor", "MovedDots", "ExpPhase" )], factor)

### Create variables needed for the accuracy analyses
# count hits, false alarms, misses, correct rejections, and RT separately for each participant (their calculation is done in Matlab: see DataProcessing.m)
data.final = ddply(data.raw,.(ParticipantNo,ExpPhase,AttendedColor,RewardedColor,MovedDots),summarize,
                  numtrials=length(which(Response!=99)), # number of trials per condition (anything that is not 99 or any other number that we're not using)
                  Hits=length(which(Response==1)), # hits: attended color moved, correct response
                  FAs=length(which(Response==2)), # false alarms: attended color did not move, (wrong) response
                  Misses=length(which(Response==0)), # misses: attended color moved, no response
                  CRs=length(which(Response==3)), # correct rejections: attended color did not move, no response
                  mean.RT=mean(RT,na.rm=TRUE)) # mean RT per condition


################################################################## Calculate accuracy and RTs per condition ###############################################################################################################################################################################################################

# Prepare the data------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

### Calculate Hits and False alarms
# Hits are calculated for each participant in each condition on trials when they are attending the color that moved. 
# False alarms are  calculated for each participant in each condition on trials when they are attending the color that didn't move (the unattended color moved, but they responded)  
# Here we create the same number of hits & fas for each of the two conditions (moved attended or not)
data.final = ddply(data.final, .(ParticipantNo,ExpPhase,AttendedColor), transform, 
                 Hits = Hits[MovedDots==AttendedColor],
                 FAs = FAs[MovedDots!=AttendedColor])

# Keep only trials on which the attended color moved (we can do behavioral analysis only on those)
data.final = subset(data.final,MovedDots==AttendedColor)

### Calculate d'
# use loglinear transformation: add 0.5 to Hits, FAs, Misses, and CRs (Hautus, 1995, Behavior Research Methods, Instruments, & Computers),
# which is preferred over the 1/2N rule (Macmillan & Kaplan, 1985, Psychological Bulletin) because it results in less biased estimates of d'.
data.final =  ddply(data.final,.(ParticipantNo,ExpPhase,RewardedColor,AttendedColor,numtrials),summarise,
                    tot.Hits=Hits+.5, # hits
                    tot.FAs=FAs+.5, # false alarms
                    tot.Misses=(numtrials-Hits)+.5, # misses
                    tot.CRs=(numtrials-FAs)+.5, # correct rejections
                    Hit.Rate=tot.Hits/(tot.Hits+tot.Misses), # hit rate
                    FA.Rate=tot.FAs/(tot.FAs+tot.CRs), # false alarm rate
                    #dprime=qnorm(Hit.Rate)-qnorm(FA.Rate),
                    Hits.RTs=mean(mean.RT,na.rm=TRUE)) # mean RTs

# Calculate SDT indices with psycho
indices = psycho::dprime(data.final$tot.Hits, data.final$tot.FAs, data.final$tot.Misses, data.final$tot.CRs) 

data.final = cbind(data.final, indices) 
                      

### Create a final dataframe for accuracy and RTs analyses
# add a new variable specifying whether the participant is attending the high or Low_Rewed color
data.final$Condition = ifelse(data.final$RewardedColor==data.final$AttendedColor,"High_Rew","Low_Rew")
# make this variable a factor for further analyses
data.final$Condition = factor(data.final$Condition)

# Exclude subjects without EEG
missing_eeg = c(5,10,13,27,14)
data.final = data.final[!data.final$ParticipantNo %in% missing_eeg,]

D prime

Means - D prime - raw data

summary = ddply(data.final,.(ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(dprime, na.rm = TRUE), digits = 2), " [", round(hdi(dprime)[[1]], digits = 2), " ", round(hdi(dprime)[[2]], digits = 2), "]")))

names(summary) = c("Reward phase", "Reward probability", "Dprime")

summary$`Reward phase` = recode(summary$`Reward phase`,
                           "Acq" = "Acquisition",
                           "Bsln" = "Baseline",
                           "Ext" = "Extinction")

summary$`Reward probability` = recode(summary$`Reward probability`,
                           "High_Rew" = "High",
                           "Low_Rew" = "Low")

summary = as.data.frame(summary)

kable(summary, caption = "Dprime per condition")
Dprime per condition
Reward phase Reward probability Dprime
Baseline High 1.64 [ -0.04 2.68 ]
Baseline Low 1.47 [ 0.04 2.3 ]
Acquisition High 1.69 [ -0.29 2.73 ]
Acquisition Low 1.62 [ 0.46 2.68 ]
Extinction High 1.6 [ -0.2 2.73 ]
Extinction Low 1.62 [ 0.74 2.88 ]

Plot

# # Plot FA rates------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Prepare the dataset
data.plot = data.final

# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] <- "Reward phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] <- "Reward probability"

# rename conditions
data.plot$`Reward phase` = recode(data.plot$`Reward phase`,
                                "Acq" = "Acquisition",
                                "Bsln" = "Baseline",
                                "Ext" = "Extinction")

data.plot$`Reward probability` = recode(data.plot$`Reward probability`,
                                      "High_Rew" = "High",
                                      "Low_Rew" = "Low")
  # Pirate plot
  pirateplot(formula=dprime ~ `Reward phase` + `Reward probability`, # dependent~independent variables
             data=data.plot, # data frame
             main="Dprime", # main title
             ylim=c(-2,4), # y-axis: limits
             ylab="Dprime", # y-axis: label
             theme=0, # preset theme (0: use your own)
             point.col="black", # points: color
             point.o=.3, # points: opacity (0-1)
             avg.line.fun = median,
             avg.line.col="black", # average line: color
             avg.line.lwd=2, # average line: line width
             avg.line.o=1, # average line: opacity (0-1)
             bean.b.col="black", # bean border, color
             bean.lwd=0.6, # bean border, line width
             bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
             bean.b.o=0.3, # bean border, opacity (0-1)
             bean.f.col="gray", # bean filling, color
             bean.f.o=.1, # bean filling, opacity (0-1)
             cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
             gl.col="gray", # gridlines: color
             gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
             cex.lab=1, # axis labels: size
             cex.axis=1, # axis numbers: size
             cex.names = 1,
             bty="l", # plot box type
             back.col="white") # background, color

Statistics

# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))


model.full.Acc.dprime = readRDS("model.full.Acc.dprime.rds")
compare.waic.Acc.dprime = readRDS("compare.Acc.dprime.waic")
compare.Acc.dprime.waic.weights = readRDS("compare.Acc.dprime.waic.weights")
bR2.null.Acc.dprime = readRDS("bR2.null.Acc.dprime")
bR2.expphase.Acc.dprime = readRDS("bR2.expphase.Acc.dprime")
bR2.full.Acc.dprime = readRDS("bR2.full.Acc.dprime")

Model comparison with WAIC

print(compare.waic.Acc.dprime)
## Output of model 'model.null.Acc.dprime':
## 
## Computed from 12000 by 258 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic   -266.6 13.2
## p_waic        30.4  3.0
## waic         533.3 26.5
## 
## Output of model 'model.expphase.Acc.dprime':
## 
## Computed from 12000 by 258 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic   -270.5 13.2
## p_waic        34.2  3.3
## waic         541.0 26.5
## 
## Output of model 'model.full.Acc.dprime':
## 
## Computed from 12000 by 258 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic   -101.3  9.9
## p_waic        77.1  4.9
## waic         202.7 19.7
## 
## Model comparisons:
##                           elpd_diff se_diff
## model.full.Acc.dprime        0.0       0.0 
## model.null.Acc.dprime     -165.3      14.5 
## model.expphase.Acc.dprime -169.2      14.5

Model comparison with WAIC weights

print(compare.Acc.dprime.waic.weights)
##     model.null.Acc.dprime model.expphase.Acc.dprime 
##              1.623994e-72              3.417889e-74 
##     model.full.Acc.dprime 
##              1.000000e+00

Bayesian R squared

Null model

print(bR2.null.Acc.dprime)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2670593 0.04908653 0.1658761 0.3580086

Experiment phase model

print(bR2.expphase.Acc.dprime)
##     Estimate  Est.Error     Q2.5     Q97.5
## R2 0.2746315 0.04876397 0.174661 0.3657918

Full model

print(bR2.full.Acc.dprime)
##    Estimate  Est.Error      Q2.5     Q97.5
## R2 0.840934 0.01411943 0.8108646 0.8657979

Checking the best model

Plotting the chains

# Plot chains
plot(model.full.Acc.dprime, pars = "^b_", ask = FALSE, N=6)

Summary of the best model

# Summary of the best model
print(tidy(model.full.Acc.dprime, par_type = "non-varying"), digits = 1)
##                           term estimate std.error lower upper
## 1                    Intercept     1.64      0.12  1.44  1.84
## 2                  ExpPhaseAcq     0.05      0.07 -0.06  0.15
## 3                  ExpPhaseExt    -0.04      0.07 -0.15  0.07
## 4             ConditionLow_Rew    -0.17      0.17 -0.45  0.11
## 5 ExpPhaseAcq:ConditionLow_Rew     0.10      0.09 -0.05  0.25
## 6 ExpPhaseExt:ConditionLow_Rew     0.18      0.09  0.04  0.33

Plotting the best model

Sample from the posterior

# Analyzing the posterior and differences between conditions

post = posterior_samples(model.full.Acc.dprime, "^b")


################################################ Baseline ####

######### High reward
Baseline_High = post[["b_Intercept"]]
######### Low reward
Baseline_Low = post[["b_Intercept"]] + 
  post[["b_ConditionLow_Rew"]] 

################################################ Acquistion

######### High reward
Acquisition_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] 
######### Low reward
Acquisition_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ExpPhaseAcq:ConditionLow_Rew"]]

################################################ Extinction

######### High reward
Extinction_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] 
######### Low reward
Extinction_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ExpPhaseExt:ConditionLow_Rew"]]
# Data from the model

# make a data frame of conditions from the posterior
posterior_conditions = melt(data.frame(Baseline_High, Baseline_Low, Acquisition_High, Acquisition_Low, Extinction_High, Extinction_Low))

posterior_conditions =  posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability"), "_", extra = "merge")

names(posterior_conditions)[3] = "dprime"

posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
                                      "High" = "High reward",
                                      "Low" = "Low reward")


posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
                                "Acquisition" = "Training",
                                "Baseline" = "Baseline",
                                "Extinction" = "Test")


# Make a table with credible intervals
posterior_means = as.data.frame(c("Baseline High Reward", 
                                  "Baseline Low Reward", 
                                  "Acquisition High Reward", 
                                  "Acquisition Low Reward", 
                                  "Extinction High Reward", 
                                  "Extinction Low Reward"))
names(posterior_means)[1] = "Condition"

posterior_means$low_95_CI = c(round(hdi(Baseline_High)[[1]], digits = 2),
                        round(hdi(Baseline_Low)[[1]], digits = 2),
                        round(hdi(Acquisition_High)[[1]], digits = 2),
                        round(hdi(Acquisition_Low)[[1]], digits = 2),
                        round(hdi(Extinction_High)[[1]], digits = 2),
                        round(hdi(Extinction_Low)[[1]], digits = 2))

posterior_means$high_95_CI = c(round(hdi(Baseline_High)[[2]], digits = 2),
                        round(hdi(Baseline_Low)[[2]], digits = 2),
                        round(hdi(Acquisition_High)[[2]], digits = 2),
                        round(hdi(Acquisition_Low)[[2]], digits = 2),
                        round(hdi(Extinction_High)[[2]], digits = 2),
                        round(hdi(Extinction_Low)[[2]], digits = 2))

posterior_means$low_50_CI = c(round(hdi(Baseline_High,0.5)[[1]], digits = 2),
                        round(hdi(Baseline_Low,0.5)[[1]], digits = 2),
                        round(hdi(Acquisition_High,0.5)[[1]], digits = 2),
                        round(hdi(Acquisition_Low,0.5)[[1]], digits = 2),
                        round(hdi(Extinction_High,0.5)[[1]], digits = 2),
                        round(hdi(Extinction_Low,0.5)[[1]], digits = 2))

posterior_means$high_50_CI = c(round(hdi(Baseline_High,0.5)[[2]], digits = 2),
                        round(hdi(Baseline_Low,0.5)[[2]], digits = 2),
                        round(hdi(Acquisition_High,0.5)[[2]], digits = 2),
                        round(hdi(Acquisition_Low,0.5)[[2]], digits = 2),
                        round(hdi(Extinction_High,0.5)[[2]], digits = 2),
                        round(hdi(Extinction_Low,0.5)[[2]], digits = 2))

posterior_means =  posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")

posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
                                      "High Reward" = "High reward",
                                      "Low Reward" = "Low reward")

posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
                                "Acquisition" = "Training",
                                "Baseline" = "Baseline",
                                "Extinction" = "Test")

# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$`dprime` = c(1,1,1,1,1,1)


# Raw data
# Prepare the dataset
data.plot = data.final

# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] = "Reward Phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] = "Reward Probability"

# rename conditions
data.plot$`Reward Phase` = recode(data.plot$`Reward Phase`,
                                "Acq" = "Training",
                                "Bsln" = "Baseline",
                                "Ext" = "Test")

data.plot$`Reward Probability` = recode(data.plot$`Reward Probability`,
                                      "High_Rew" = "High reward",
                                      "Low_Rew" = "Low reward")


# Violin + the model
# tiff("dprime.tiff", units="in", width=8, height=5, res=800)

f.2.a = ggplot(data=data.plot, aes(x=`Reward Phase`, y=`dprime`, label=`Reward Probability`)) + 
  # violin of the raw data
  geom_violin(data=data.plot, color ="#636363",trim = TRUE,alpha="0.2") +
  # individual participants for the raw data
  # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
  geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
  # mean of the model
  stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`dprime`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
  # 95 credible intervals of the model
  geom_linerange(data = posterior_means,width=.1, aes(ymin=`dprime`-`dprime`+low_95_CI, ymax=`dprime`-`dprime`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
  # 50 credible intervals of the model
  geom_linerange(data = posterior_means,width=.1, aes(ymin=`dprime`-`dprime`+low_50_CI, ymax=`dprime`-`dprime`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
  scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
  facet_grid(.~`Reward Probability`) + 
  scale_y_continuous(limits = c(-1,3),breaks=seq(-1,3,0.5)) +
  labs(title="Sensitivity d′", y = "Sensitivity d′") + 
  theme(
    axis.line = element_line(size=1, colour = "black"),
    panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
    panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
    panel.border = element_blank(), 
    panel.background = element_blank(),
    plot.title = element_text(size = 18,  face = "bold",hjust = 0.5),
    text=element_text(colour="black", size = 14),
    axis.text.x=element_text(colour="black", size = 14),
    axis.text.y=element_text(colour="black", size = 14),
    axis.title=element_text(size=16,colour = "black")) + 
  theme(strip.background =element_rect(fill="white")) + 
  theme(strip.text = element_text(size = 16))
f.2.a

# dev.off()

Table of means across conditions

library(HDInterval)

# Make a table with conditions
posterior_means = as.data.frame(c("Baseline High Reward", 
                                  "Baseline Low Reward", 
                                  "Acquisition High Reward", 
                                  "Acquisition Low Reward", 
                                  "Extinction High Reward", 
                                  "Extinction Low Reward"))
names(posterior_means)[1] = "Condition"

posterior_means$Mean = c(paste(round(mean(Baseline_High), digits = 2), " [", round(hdi(Baseline_High)[[1]], digits = 2), " ", round(hdi(Baseline_High)[[2]], digits = 2), "]"),
                        paste(round(mean(Baseline_Low), digits = 2), " [", round(hdi(Baseline_Low)[[1]], digits = 2), " ", round(hdi(Baseline_Low)[[2]], digits = 2), "]"),
                        paste(round(mean(Acquisition_High), digits = 2), " [", round(hdi(Acquisition_High)[[1]], digits = 2), " ", round(hdi(Acquisition_High)[[2]], digits = 2), "]"),
                        paste(round(mean(Acquisition_Low), digits = 2), " [", round(hdi(Acquisition_Low)[[1]], digits = 2), " ", round(hdi(Acquisition_Low)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_High), digits = 2), " [", round(hdi(Extinction_High)[[1]], digits = 2), " ", round(hdi(Extinction_High)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_Low), digits = 2), " [", round(hdi(Extinction_Low)[[1]], digits = 2), " ", round(hdi(Extinction_Low)[[2]], digits = 2), "]"))

names(posterior_means)[2] = "Mean [HDI]"

posterior_means =  posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")
                        
kable(posterior_means, caption = "Means per condition")
Means per condition
Reward Phase Reward Probability Mean [HDI]
Baseline High Reward 1.64 [ 1.39 1.87 ]
Baseline Low Reward 1.48 [ 1.25 1.69 ]
Acquisition High Reward 1.69 [ 1.44 1.93 ]
Acquisition Low Reward 1.62 [ 1.41 1.84 ]
Extinction High Reward 1.61 [ 1.36 1.84 ]
Extinction Low Reward 1.62 [ 1.41 1.84 ]

Inference about the best model

Check the difference between high reward in baseline vs. acquisition

Diff_Bsln_Acq_High = Acquisition_High - Baseline_High
plotPost(Diff_Bsln_Acq_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

length(Diff_Bsln_Acq_High)
## [1] 12000

Check the difference between low reward in baseline vs. acquisition

Diff_Bsln_Acq_Low = Acquisition_Low - Baseline_Low
plotPost(Diff_Bsln_Acq_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

length(Diff_Bsln_Acq_Low)
## [1] 12000

Check the difference between high and low reward in baseline vs. acquisition

Diff_Bsln_Acq_High_vs_Low = Diff_Bsln_Acq_Low - Diff_Bsln_Acq_High 
plotPost(Diff_Bsln_Acq_High_vs_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

paste("Mean = ",round(mean(Diff_Bsln_Acq_High_vs_Low), digits = 2), " [", round(hdi(Diff_Bsln_Acq_High_vs_Low)[[1]], digits = 2), " ", round(hdi(Diff_Bsln_Acq_High_vs_Low)[[2]], digits = 2), "]")
## [1] "Mean =  0.1  [ -0.08   0.27 ]"
length(Diff_Bsln_Acq_High_vs_Low)
## [1] 12000

Check the difference between high reward in acquisition vs. extinction

Diff_Acq_Ext_High = Extinction_High - Acquisition_High
plotPost(Diff_Acq_Ext_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

length(Diff_Acq_Ext_High)
## [1] 12000

Check the difference between low reward in acquisition vs. extinction

Diff_Acq_Ext_Low = Extinction_Low - Acquisition_Low
plotPost(Diff_Acq_Ext_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

length(Diff_Acq_Ext_Low)
## [1] 12000

Reaction times

Means - raw data

summary = ddply(data.final,.(ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(Hits.RTs, na.rm = TRUE), digits = 2), " [", round(hdi(Hits.RTs)[[1]], digits = 2), " ", round(hdi(Hits.RTs)[[2]], digits = 2), "]")))

names(summary) = c("Reward phase", "Reward probability", "Reaction times")

summary$`Reward phase` = recode(summary$`Reward phase`,
                           "Acq" = "Acquisition",
                           "Bsln" = "Baseline",
                           "Ext" = "Extinction")

summary$`Reward probability` = recode(summary$`Reward probability`,
                           "High_Rew" = "High",
                           "Low_Rew" = "Low")

summary = as.data.frame(summary)

kable(summary, caption = "Reaction times per condition")
Reaction times per condition
Reward phase Reward probability Reaction times
Baseline High 546.59 [ 485.64 619.34 ]
Baseline Low 551.1 [ 490.5 631.36 ]
Acquisition High 524.99 [ 467.12 599.49 ]
Acquisition Low 537.94 [ 465.32 584.63 ]
Extinction High 528.98 [ 457.08 599.83 ]
Extinction Low 539.75 [ 455.8 623.21 ]

Plot

# Prepare the dataset
data.plot = data.final

# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] = "Reward phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] = "Reward probability"

# rename conditions
data.plot$`Reward phase` = recode(data.plot$`Reward phase`,
                                "Acq" = "Acquisition",
                                "Bsln" = "Baseline",
                                "Ext" = "Extinction")

data.plot$`Reward probability` = recode(data.plot$`Reward probability`,
                                      "High_Rew" = "High",
                                      "Low_Rew" = "Low")

# Pirate plot
  pirateplot(formula = Hits.RTs ~ `Reward phase` + `Reward probability`, # dependent~independent variables
             data=data.plot, # data frame
             main="Reaction times", # main title
             ylim=c(400,700), # y-axis: limits
             ylab="Reaction time", # y-axis: label
             theme=0, # preset theme (0: use your own)
             point.col="black", # points: color
             point.o=.3, # points: opacity (0-1)
             avg.line.fun = median,
             avg.line.col="black", # average line: color
             avg.line.lwd=2, # average line: line width
             avg.line.o=1, # average line: opacity (0-1)
             bean.b.col="black", # bean border, color
             bean.lwd=0.6, # bean border, line width
             bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
             bean.b.o=0.3, # bean border, opacity (0-1)
             bean.f.col="gray", # bean filling, color
             bean.f.o=.1, # bean filling, opacity (0-1)
             cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
             gl.col="gray", # gridlines: color
             gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
             cex.lab=1, # axis labels: size
             cex.axis=1, # axis numbers: size
             cex.names = 1,
             bty="l", # plot box type
             back.col="white") # background, color

Statistics

# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))

model.full.RT = readRDS("model.full.RT.rds")
compare.waic.RT = readRDS("compare.RT.waic")
compare.RT.waic.weights = readRDS("compare.RT.waic.weights")
bR2.null.RT = readRDS("bR2.null.RT")
bR2.expphase.RT = readRDS("bR2.expphase.RT")
bR2.full.RT = readRDS("bR2.full.RT")

Model comparison with WAIC

print(compare.waic.RT)
## Output of model 'model.null.RT':
## 
## Computed from 12000 by 258 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic  -1250.1 15.8
## p_waic        36.4  4.6
## waic        2500.2 31.6
## 
## Output of model 'model.expphase.RT':
## 
## Computed from 12000 by 258 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic  -1241.5 17.6
## p_waic        48.7  7.3
## waic        2483.0 35.3
## 
## Output of model 'model.full.RT':
## 
## Computed from 12000 by 258 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic  -1161.3 15.0
## p_waic        93.4  9.4
## waic        2322.5 30.0
## 
## Model comparisons:
##                   elpd_diff se_diff
## model.full.RT       0.0       0.0  
## model.expphase.RT -80.2       8.9  
## model.null.RT     -88.8       9.0

Bayesian R squared

Null model

print(bR2.null.RT)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.5016012 0.03697726 0.4229632 0.5683283

Experiment phase model

print(bR2.expphase.RT)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.5655434 0.03698904 0.4884728 0.6331867

Full model

print(bR2.full.RT)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.8168505 0.02262293 0.7666812 0.8545284

Checking the best model

Plotting the chains

# Plot chains
plot(model.full.RT, pars = "^b_", ask = FALSE, N=6)

Summary of the best model

# Summary of the best model
print(tidy(model.full.RT, par_type = "non-varying"), digits = 1)
##                           term estimate std.error lower upper
## 1                    Intercept      547         6 536.2   557
## 2                  ExpPhaseAcq      -22         4 -28.8   -14
## 3                  ExpPhaseExt      -18         6 -26.6    -8
## 4             ConditionLow_Rew        5         6  -5.4    15
## 5 ExpPhaseAcq:ConditionLow_Rew        8         5  -0.3    17
## 6 ExpPhaseExt:ConditionLow_Rew        6         5  -2.4    15

Plotting the best model

Sample from the posterior

# Analyzing the posterior and differences between conditions

post = posterior_samples(model.full.RT, "^b")


################################################ Baseline ####

######### High reward
Baseline_High = post[["b_Intercept"]]
######### Low reward
Baseline_Low = post[["b_Intercept"]] + 
  post[["b_ConditionLow_Rew"]] 

################################################ Acquistion

######### High reward
Acquisition_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] 
######### Low reward
Acquisition_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ExpPhaseAcq:ConditionLow_Rew"]]

################################################ Extinction

######### High reward
Extinction_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] 
######### Low reward
Extinction_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ExpPhaseExt:ConditionLow_Rew"]]
# Data from the model

# make a data frame of conditions from the posterior
posterior_conditions = melt(data.frame(Baseline_High, Baseline_Low, Acquisition_High, Acquisition_Low, Extinction_High, Extinction_Low))

posterior_conditions =  posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability"), "_", extra = "merge")

names(posterior_conditions)[3] = "Reaction time"

posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
                                      "High" = "High reward",
                                      "Low" = "Low reward")

posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
                                "Acquisition" = "Training",
                                "Baseline" = "Baseline",
                                "Extinction" = "Test")


# Make a table with credible intervals
posterior_means = as.data.frame(c("Baseline High Reward", 
                                  "Baseline Low Reward", 
                                  "Acquisition High Reward", 
                                  "Acquisition Low Reward", 
                                  "Extinction High Reward", 
                                  "Extinction Low Reward"))
names(posterior_means)[1] = "Condition"

posterior_means$low_95_CI = c(round(hdi(Baseline_High)[[1]], digits = 2),
                        round(hdi(Baseline_Low)[[1]], digits = 2),
                        round(hdi(Acquisition_High)[[1]], digits = 2),
                        round(hdi(Acquisition_Low)[[1]], digits = 2),
                        round(hdi(Extinction_High)[[1]], digits = 2),
                        round(hdi(Extinction_Low)[[1]], digits = 2))

posterior_means$high_95_CI = c(round(hdi(Baseline_High)[[2]], digits = 2),
                        round(hdi(Baseline_Low)[[2]], digits = 2),
                        round(hdi(Acquisition_High)[[2]], digits = 2),
                        round(hdi(Acquisition_Low)[[2]], digits = 2),
                        round(hdi(Extinction_High)[[2]], digits = 2),
                        round(hdi(Extinction_Low)[[2]], digits = 2))

posterior_means$low_50_CI = c(round(hdi(Baseline_High,0.5)[[1]], digits = 2),
                        round(hdi(Baseline_Low,0.5)[[1]], digits = 2),
                        round(hdi(Acquisition_High,0.5)[[1]], digits = 2),
                        round(hdi(Acquisition_Low,0.5)[[1]], digits = 2),
                        round(hdi(Extinction_High,0.5)[[1]], digits = 2),
                        round(hdi(Extinction_Low,0.5)[[1]], digits = 2))

posterior_means$high_50_CI = c(round(hdi(Baseline_High,0.5)[[2]], digits = 2),
                        round(hdi(Baseline_Low,0.5)[[2]], digits = 2),
                        round(hdi(Acquisition_High,0.5)[[2]], digits = 2),
                        round(hdi(Acquisition_Low,0.5)[[2]], digits = 2),
                        round(hdi(Extinction_High,0.5)[[2]], digits = 2),
                        round(hdi(Extinction_Low,0.5)[[2]], digits = 2))

posterior_means =  posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")

posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
                                      "High Reward" = "High reward",
                                      "Low Reward" = "Low reward")

posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
                                "Acquisition" = "Training",
                                "Baseline" = "Baseline",
                                "Extinction" = "Test")

# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$`Reaction time` = c(555,555,555,555,555,555)


# Raw data
# Prepare the dataset
data.plot = data.final

# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] = "Reward Phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] = "Reward Probability"

# rename conditions
data.plot$`Reward Phase` = recode(data.plot$`Reward Phase`,
                                "Acq" = "Training",
                                "Bsln" = "Baseline",
                                "Ext" = "Test")

data.plot$`Reward Probability` = recode(data.plot$`Reward Probability`,
                                      "High_Rew" = "High reward",
                                      "Low_Rew" = "Low reward")

data.plot$`Reaction time` = data.plot$Hits.RTs


# Violin + the model
# tiff("RTs.tiff", units="in", width=8, height=5, res=800)

f.2.b = ggplot(data=data.plot, aes(x=`Reward Phase`, y=`Reaction time`, label=`Reward Probability`)) + 
  # violin of the raw data
  geom_violin(data=data.plot, color ="#636363",trim = TRUE,alpha="0.2") +
  # individual participants for the raw data
  # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
  geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
  # mean of the model
  stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`Reaction time`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
  # 95 credible intervals of the model
  geom_linerange(data = posterior_means,width=.1, aes(ymin=`Reaction time`-`Reaction time`+low_95_CI, ymax=`Reaction time`-`Reaction time`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
  # 50 credible intervals of the model
  geom_linerange(data = posterior_means,width=.1, aes(ymin=`Reaction time`-`Reaction time`+low_50_CI, ymax=`Reaction time`-`Reaction time`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
  scale_x_discrete(limits=c("Baseline", "Training", "Test")) +
  facet_grid(.~`Reward Probability`) + 
  scale_y_continuous(limits = c(450,650), breaks=seq(450,650,50)) +
  labs(title="Reaction time", y = "Reaction time (ms)") + 
  theme(
    axis.line = element_line(size=1, colour = "black"),
    panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
    panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
    panel.border = element_blank(), 
    panel.background = element_blank(),
    plot.title = element_text(size = 18,  face = "bold",hjust = 0.5),
    text=element_text(colour="black", size = 14),
    axis.text.x=element_text(colour="black", size = 14),
    axis.text.y=element_text(colour="black", size = 14),
    axis.title=element_text(size=16,colour = "black")) + 
  theme(strip.background =element_rect(fill="white")) + 
  theme(strip.text = element_text(size = 16))
f.2.b 

# dev.off()

Table of means across conditions

library(HDInterval)

# Make a table with conditions
posterior_means = as.data.frame(c("Baseline High Reward", 
                                  "Baseline Low Reward", 
                                  "Acquisition High Reward", 
                                  "Acquisition Low Reward", 
                                  "Extinction High Reward", 
                                  "Extinction Low Reward"))
names(posterior_means)[1] = "Condition"

posterior_means$Mean = c(paste(round(mean(Baseline_High), digits = 2), " [", round(hdi(Baseline_High)[[1]], digits = 2), " ", round(hdi(Baseline_High)[[2]], digits = 2), "]"),
                        paste(round(mean(Baseline_Low), digits = 2), " [", round(hdi(Baseline_Low)[[1]], digits = 2), " ", round(hdi(Baseline_Low)[[2]], digits = 2), "]"),
                        paste(round(mean(Acquisition_High), digits = 2), " [", round(hdi(Acquisition_High)[[1]], digits = 2), " ", round(hdi(Acquisition_High)[[2]], digits = 2), "]"),
                        paste(round(mean(Acquisition_Low), digits = 2), " [", round(hdi(Acquisition_Low)[[1]], digits = 2), " ", round(hdi(Acquisition_Low)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_High), digits = 2), " [", round(hdi(Extinction_High)[[1]], digits = 2), " ", round(hdi(Extinction_High)[[2]], digits = 2), "]"),
                        paste(round(mean(Extinction_Low), digits = 2), " [", round(hdi(Extinction_Low)[[1]], digits = 2), " ", round(hdi(Extinction_Low)[[2]], digits = 2), "]"))

names(posterior_means)[2] = "Mean [HDI]"

posterior_means =  posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")

kable(posterior_means, caption = "Means per condition")
Means per condition
Reward Phase Reward Probability Mean [HDI]
Baseline High Reward 546.54 [ 534.33 559.3 ]
Baseline Low Reward 551.13 [ 539.34 563.5 ]
Acquisition High Reward 524.91 [ 512.94 536.3 ]
Acquisition Low Reward 537.99 [ 526.48 550.32 ]
Extinction High Reward 528.97 [ 515.9 541.99 ]
Extinction Low Reward 539.85 [ 525.63 554.34 ]

Inference about the best model

Check the difference between high reward in baseline vs. acquisition

Diff_Bsln_Acq_High = Acquisition_High - Baseline_High
plotPost(Diff_Bsln_Acq_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

length(Diff_Bsln_Acq_High)
## [1] 12000

Check the difference between low reward in baseline vs. acquisition

Diff_Bsln_Acq_Low = Acquisition_Low - Baseline_Low 
plotPost(Diff_Bsln_Acq_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

length(Diff_Bsln_Acq_Low)
## [1] 12000

Check the difference between high and low reward in baseline vs. acquisition

Diff_Bsln_Acq_High_vs_Low = Diff_Bsln_Acq_High - Diff_Bsln_Acq_Low 
plotPost(Diff_Bsln_Acq_High_vs_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

paste("Mean = ",round(mean(Diff_Bsln_Acq_High_vs_Low), digits = 2), " [", round(hdi(Diff_Bsln_Acq_High_vs_Low)[[1]], digits = 2), " ", round(hdi(Diff_Bsln_Acq_High_vs_Low)[[2]], digits = 2), "]")
## [1] "Mean =  -8.49  [ -18.61   2.06 ]"
length(Diff_Bsln_Acq_High_vs_Low)
## [1] 12000

Check the difference between high reward in acquisition vs. extinction

Diff_Acq_Ext_High = Extinction_High - Acquisition_High
plotPost(Diff_Acq_Ext_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

length(Diff_Acq_Ext_High)
## [1] 12000

Check the difference between low reward in acquisition vs. extinction

Diff_Acq_Ext_Low = Extinction_Low - Acquisition_Low
plotPost(Diff_Acq_Ext_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

length(Diff_Acq_Ext_Low)
## [1] 12000

Behavior - splitting the phases

Import data

# Clear environemnt and import data------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

# # clear the environment
# rm(list=ls()) 
# # clear the console
# cat("\014") 
# #load packages and install them if they're not installed
# if (!require("pacman")) install.packages("pacman")
# pacman::p_load(plyr,Rmisc,yarrr,BayesFactor,reshape2,brms, broom, tidyverse, brmstools, BEST, knitr, here,psych)
# # set seed
# set.seed(42) 
# # set directory
# setwd(here())
# import data
data.raw = read.csv(file = here("data","Data_behavior_exp1_48pps.csv"),header=TRUE,na.strings="NaN")

# Prepare the dataset------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

### Adding and renaming variables 
# rename EventType variable
names(data.raw)[names(data.raw) == "EventType"] = "MovedDots" 
# add a variable with the name of the attended color instead of a numbers
data.raw$AttendedColor = ifelse(data.raw$AttendedColor==1,"red","blue")
# add a variable saying which color was linked with High_Rew (even numbers - blue was High_Rew)
data.raw$RewardedColor = ifelse(data.raw$ParticipantNo%%2==0,"blue","red") 
# add a variable with the name of the moved color instead of a numbers
data.raw$MovedDots = ifelse(data.raw$MovedDots==1,"red","blue") 
# split experimental phases into 6 isntead of 3 phases (trial 0-200: Bsln; trial 201-400: Acq; trial 401-600: Ext)
data.raw$ExpPhase = cut(data.raw$Trial,breaks=c(0,100,200,300,400,500,600),labels=c("Bsln1","Bsln2","Acq1","Acq2","Ext1","Ext2"))
# split experimental phases into 3 phases (trial 0-200: Bsln; trial 201-400: Acq; trial 401-600: Ext)
# data.raw$ExpPhase = cut(data.raw$Trial,breaks=c(0,200,400,600),labels=c("Bsln","Acq","Ext")) # trial 0-200: Bsln; trial 201-400: Acq; trial 401-600: Ext

### Convert variables to be used in analyses into factors
data.raw[c("ParticipantNo", "AttendedColor","RewardedColor", "MovedDots", "ExpPhase" )] = 
  lapply(data.raw[c("ParticipantNo", "AttendedColor","RewardedColor", "MovedDots", "ExpPhase" )], factor)

### Create variables needed for the accuracy analyses
# count hits, false alarms, misses, correct rejections, and RT separately for each participant (their calculation is done in Matlab: see DataProcessing.m)
data.final = ddply(data.raw,.(ParticipantNo,ExpPhase,AttendedColor,RewardedColor,MovedDots),summarize,
                  numtrials=length(which(Response!=99)), # number of trials per condition (anything that is not 99 or any other number that we're not using)
                  Hits=length(which(Response==1)), # hits: attended color moved, correct response
                  FAs=length(which(Response==2)), # false alarms: attended color did not move, (wrong) response
                  Misses=length(which(Response==0)), # misses: attended color moved, no response
                  CRs=length(which(Response==3)), # correct rejections: attended color did not move, no response
                  mean.RT=mean(RT,na.rm=TRUE)) # mean RT per condition


################################################################## Calculate accuracy and RTs per condition ###############################################################################################################################################################################################################

# Prepare the data------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------

### Calculate Hits and False alarms
# Hits are calculated for each participant in each condition on trials when they are attending the color that moved. 
# False alarms are  calculated for each participant in each condition on trials when they are attending the color that didn't move (the unattended color moved, but they responded)  
# Here we create the same number of hits & fas for each of the two conditions (moved attended or not)
data.final = ddply(data.final, .(ParticipantNo,ExpPhase,AttendedColor), transform, 
                 Hits = Hits[MovedDots==AttendedColor],
                 FAs = FAs[MovedDots!=AttendedColor])

# Keep only trials on which the attended color moved (we can do behavioral analysis only on those)
data.final = subset(data.final,MovedDots==AttendedColor)

### Calculate d'
# use loglinear transformation: add 0.5 to Hits, FAs, Misses, and CRs (Hautus, 1995, Behavior Research Methods, Instruments, & Computers),
# which is preferred over the 1/2N rule (Macmillan & Kaplan, 1985, Psychological Bulletin) because it results in less biased estimates of d'.
data.final =  ddply(data.final,.(ParticipantNo,ExpPhase,RewardedColor,AttendedColor,numtrials),summarise,
                    tot.Hits=Hits+.5, # hits
                    tot.FAs=FAs+.5, # false alarms
                    tot.Misses=(numtrials-Hits)+.5, # misses
                    tot.CRs=(numtrials-FAs)+.5, # correct rejections
                    Hit.Rate=tot.Hits/(tot.Hits+tot.Misses), # hit rate
                    FA.Rate=tot.FAs/(tot.FAs+tot.CRs), # false alarm rate
                    #dprime=qnorm(Hit.Rate)-qnorm(FA.Rate),
                    Hits.RTs=mean(mean.RT,na.rm=TRUE)) # mean RTs

# Calculate SDT indices with psycho
indices = psycho::dprime(data.final$tot.Hits, data.final$tot.FAs, data.final$tot.Misses, data.final$tot.CRs) 

data.final = cbind(data.final, indices) 
                      

### Create a final dataframe for accuracy and RTs analyses
# add a new variable specifying whether the participant is attending the high or Low_Rewed color
data.final$Condition = ifelse(data.final$RewardedColor==data.final$AttendedColor,"High_Rew","Low_Rew")
# make this variable a factor for further analyses
data.final$Condition = factor(data.final$Condition)

# Exclude subjects without EEG
missing_eeg = c(5,10,13,27,14)
data.final = data.final[!data.final$ParticipantNo %in% missing_eeg,]

D prime

Means - D prime - raw data

summary = ddply(data.final,.(ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(dprime, na.rm = TRUE), digits = 2), " [", round(hdi(dprime)[[1]], digits = 2), " ", round(hdi(dprime)[[2]], digits = 2), "]")))

names(summary) = c("Reward phase", "Reward probability", "Dprime")

# rename conditions
summary$`Reward phase` = recode(summary$`Reward phase`,
                                "Acq1" = "Acquisition1",
                                "Acq2" = "Acquisition2",
                                "Bsln1" = "Baseline1",
                                "Bsln2" = "Baseline2",
                                "Ext1" = "Extinction1",
                                "Ext2" = "Extinction2")

summary$`Reward probability` = recode(summary$`Reward probability`,
                           "High_Rew" = "High",
                           "Low_Rew" = "Low")

summary = as.data.frame(summary)

kable(summary, caption = "Dprime per condition")
Dprime per condition
Reward phase Reward probability Dprime
Baseline1 High 1.48 [ -0.36 2.62 ]
Baseline1 Low 1.32 [ 0.09 2.35 ]
Baseline2 High 1.6 [ -0.27 2.56 ]
Baseline2 Low 1.47 [ 0.08 2.33 ]
Acquisition1 High 1.54 [ -0.08 2.65 ]
Acquisition1 Low 1.59 [ 0.47 2.45 ]
Acquisition2 High 1.59 [ 0.08 2.56 ]
Acquisition2 Low 1.48 [ 0 2.62 ]
Extinction1 High 1.48 [ -0.07 2.47 ]
Extinction1 Low 1.5 [ 0.36 2.5 ]
Extinction2 High 1.49 [ -0.38 2.49 ]
Extinction2 Low 1.55 [ 0.65 2.55 ]

Plot

# # Plot FA rates------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------
# Prepare the dataset
data.plot = data.final

# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] <- "Reward phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] <- "Reward probability"

# rename conditions
data.plot$`Reward phase` = recode(data.plot$`Reward phase`,
                                "Acq1" = "Acquisition1",
                                "Acq2" = "Acquisition2",
                                "Bsln1" = "Baseline1",
                                "Bsln2" = "Baseline2",
                                "Ext1" = "Extinction1",
                                "Ext2" = "Extinction2")

data.plot$`Reward probability` = recode(data.plot$`Reward probability`,
                                      "High_Rew" = "High",
                                      "Low_Rew" = "Low")
  # Pirate plot
  pirateplot(formula=dprime ~ `Reward phase` + `Reward probability`, # dependent~independent variables
             data=data.plot, # data frame
             main="Dprime", # main title
             ylim=c(-2,4), # y-axis: limits
             ylab="Dprime", # y-axis: label
             theme=0, # preset theme (0: use your own)
             point.col="black", # points: color
             point.o=.3, # points: opacity (0-1)
             avg.line.fun = median,
             avg.line.col="black", # average line: color
             avg.line.lwd=2, # average line: line width
             avg.line.o=1, # average line: opacity (0-1)
             bean.b.col="black", # bean border, color
             bean.lwd=0.6, # bean border, line width
             bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
             bean.b.o=0.3, # bean border, opacity (0-1)
             bean.f.col="gray", # bean filling, color
             bean.f.o=.1, # bean filling, opacity (0-1)
             cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
             gl.col="gray", # gridlines: color
             gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
             cex.lab=1, # axis labels: size
             cex.axis=1, # axis numbers: size
             cex.names = 1,
             bty="l", # plot box type
             back.col="white") # background, color

Statistics

# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))


model.full.Acc.dprime.split = readRDS("model.full.Acc.dprime.split.rds")
compare.waic.Acc.dprime.split = readRDS("compare.Acc.dprime.split.waic")
compare.Acc.dprime.split.waic.weights = readRDS("compare.Acc.dprime.split.waic.weights")
bR2.null.Acc.dprime.split = readRDS("bR2.null.Acc.dprime.split")
bR2.expphase.Acc.dprime.split = readRDS("bR2.expphase.Acc.dprime.split")
bR2.full.Acc.dprime.split = readRDS("bR2.full.Acc.dprime.split")

Model comparison with WAIC

print(compare.waic.Acc.dprime.split)
## Output of model 'model.null.Acc.dprime':
## 
## Computed from 12000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic   -500.8 17.7
## p_waic        36.2  2.5
## waic        1001.5 35.4
## 
## Output of model 'model.expphase.Acc.dprime':
## 
## Computed from 12000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic   -509.5 17.7
## p_waic        48.1  3.2
## waic        1019.1 35.3
## 
## Output of model 'model.full.Acc.dprime':
## 
## Computed from 12000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic   -275.4 17.0
## p_waic       102.8  6.3
## waic         550.8 34.1
## 
## Model comparisons:
##                           elpd_diff se_diff
## model.full.Acc.dprime        0.0       0.0 
## model.null.Acc.dprime     -225.4      19.9 
## model.expphase.Acc.dprime -234.1      19.9

Model comparison with WAIC weights

print(compare.Acc.dprime.split.waic.weights)
##     model.null.Acc.dprime model.expphase.Acc.dprime 
##              1.343939e-98             2.079420e-102 
##     model.full.Acc.dprime 
##              1.000000e+00

Bayesian R squared

Null model

print(bR2.null.Acc.dprime.split)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.2864918 0.03099612 0.2239089 0.3449991

Experiment phase model

print(bR2.expphase.Acc.dprime.split)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.3007569 0.03047797 0.2382382 0.3576129

Full model

print(bR2.full.Acc.dprime.split)
##     Estimate  Est.Error      Q2.5    Q97.5
## R2 0.7493958 0.01352598 0.7212571 0.774454

Checking the best model

Plotting the chains

# Plot chains
plot(model.full.Acc.dprime.split, pars = "^b_", ask = FALSE, N=6)

Summary of the best model

# Summary of the best model
print(tidy(model.full.Acc.dprime.split, par_type = "non-varying"), digits = 1)
##                              term estimate std.error  lower upper
## 1                       Intercept    1.484      0.12  1.289   1.7
## 2                   ExpPhaseBsln2    0.116      0.08 -0.023   0.3
## 3                    ExpPhaseAcq1    0.060      0.08 -0.074   0.2
## 4                    ExpPhaseAcq2    0.113      0.08 -0.019   0.2
## 5                    ExpPhaseExt1    0.004      0.08 -0.128   0.1
## 6                    ExpPhaseExt2    0.011      0.08 -0.124   0.1
## 7                ConditionLow_Rew   -0.164      0.17 -0.439   0.1
## 8  ExpPhaseBsln2:ConditionLow_Rew    0.040      0.11 -0.145   0.2
## 9   ExpPhaseAcq1:ConditionLow_Rew    0.208      0.11  0.024   0.4
## 10  ExpPhaseAcq2:ConditionLow_Rew    0.048      0.11 -0.136   0.2
## 11  ExpPhaseExt1:ConditionLow_Rew    0.175      0.11 -0.009   0.4
## 12  ExpPhaseExt2:ConditionLow_Rew    0.215      0.11  0.033   0.4

Plotting the best model

Sample from the posterior

# Analyzing the posterior and differences between conditions

post = posterior_samples(model.full.Acc.dprime.split, "^b")


################################################ Baseline 1 ####

######### High reward
Baseline1_High = post[["b_Intercept"]]
######### Low reward
Baseline1_Low = post[["b_Intercept"]] + 
  post[["b_ConditionLow_Rew"]] 

################################################ Baseline 2 ####

######### High reward
Baseline2_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseBsln2"]] 
######### Low reward
Baseline2_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseBsln2"]] +
  post[["b_ConditionLow_Rew"]] +
  post[["b_ExpPhaseBsln2:ConditionLow_Rew"]]

################################################ Acquistion 1

######### High reward
Acquisition1_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq1"]] 
######### Low reward
Acquisition1_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq1"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ExpPhaseAcq1:ConditionLow_Rew"]]

################################################ Acquistion 2

######### High reward
Acquisition2_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq2"]] 
######### Low reward
Acquisition2_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq2"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ExpPhaseAcq2:ConditionLow_Rew"]]

################################################ Extinction 1

######### High reward
Extinction1_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt1"]] 
######### Low reward
Extinction1_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt1"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ExpPhaseExt1:ConditionLow_Rew"]]

################################################ Extinction 2

######### High reward
Extinction2_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt2"]] 
######### Low reward
Extinction2_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt2"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ExpPhaseExt2:ConditionLow_Rew"]]
# Data from the model

# make a data frame
posterior_conditions = melt(data.frame(Baseline1_High, Baseline2_High, Baseline1_Low, Baseline2_Low, Acquisition1_High, Acquisition1_Low, Acquisition2_High, Acquisition2_Low, Extinction1_High, Extinction1_Low, Extinction2_High, Extinction2_Low))

posterior_conditions =  posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability"), "_", extra = "merge")

names(posterior_conditions)[3] = "dprime"

posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
                                      "High" = "High reward",
                                      "Low" = "Low reward")

posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
                                "Acquisition1" = "Training1",
                                "Acquisition2" = "Training2",
                                "Baseline1" = "Baseline1",
                                "Baseline1" = "Baseline2",
                                "Extinction1" = "Test1",
                                "Extinction2" = "Test2")

# Make a table with credible intervals
posterior_means = as.data.frame(c("Baseline1 High Reward", 
                                  "Baseline2 High Reward", 
                                  "Baseline1 Low Reward", 
                                  "Baseline2 Low Reward",
                                  "Acquisition1 High Reward", 
                                  "Acquisition2 High Reward",
                                  "Acquisition1 Low Reward", 
                                  "Acquisition2 Low Reward", 
                                  "Extinction1 High Reward", 
                                  "Extinction2 High Reward",
                                  "Extinction1 Low Reward",
                                  "Extinction2 Low Reward"))
names(posterior_means)[1] = "Condition"

posterior_means$low_95_CI = c(round(hdi(Baseline1_High)[[1]], digits = 2),
                              round(hdi(Baseline2_High)[[1]], digits = 2),
                        round(hdi(Baseline1_Low)[[1]], digits = 2),
                        round(hdi(Baseline2_Low)[[1]], digits = 2),
                        round(hdi(Acquisition1_High)[[1]], digits = 2),
                        round(hdi(Acquisition2_High)[[1]], digits = 2),
                        round(hdi(Acquisition1_Low)[[1]], digits = 2),
                        round(hdi(Acquisition2_Low)[[1]], digits = 2),
                        round(hdi(Extinction1_High)[[1]], digits = 2),
                        round(hdi(Extinction2_High)[[1]], digits = 2),
                        round(hdi(Extinction1_Low)[[1]], digits = 2),
                        round(hdi(Extinction2_Low)[[1]], digits = 2))

posterior_means$high_95_CI = c(round(hdi(Baseline1_High)[[2]], digits = 2),
                              round(hdi(Baseline2_High)[[2]], digits = 2),
                        round(hdi(Baseline1_Low)[[2]], digits = 2),
                        round(hdi(Baseline2_Low)[[2]], digits = 2),
                        round(hdi(Acquisition1_High)[[2]], digits = 2),
                        round(hdi(Acquisition2_High)[[2]], digits = 2),
                        round(hdi(Acquisition1_Low)[[2]], digits = 2),
                        round(hdi(Acquisition2_Low)[[2]], digits = 2),
                        round(hdi(Extinction1_High)[[2]], digits = 2),
                        round(hdi(Extinction2_High)[[2]], digits = 2),
                        round(hdi(Extinction1_Low)[[2]], digits = 2),
                        round(hdi(Extinction2_Low)[[2]], digits = 2))

posterior_means$low_50_CI = c(round(hdi(Baseline1_High,0.5)[[1]], digits = 2),
                              round(hdi(Baseline2_High,0.5)[[1]], digits = 2),
                        round(hdi(Baseline1_Low,0.5)[[1]], digits = 2),
                        round(hdi(Baseline2_Low,0.5)[[1]], digits = 2),
                        round(hdi(Acquisition1_High,0.5)[[1]], digits = 2),
                        round(hdi(Acquisition2_High,0.5)[[1]], digits = 2),
                        round(hdi(Acquisition1_Low,0.5)[[1]], digits = 2),
                        round(hdi(Acquisition2_Low,0.5)[[1]], digits = 2),
                        round(hdi(Extinction1_High,0.5)[[1]], digits = 2),
                        round(hdi(Extinction2_High,0.5)[[1]], digits = 2),
                        round(hdi(Extinction1_Low,0.5)[[1]], digits = 2),
                        round(hdi(Extinction2_Low,0.5)[[1]], digits = 2))

posterior_means$high_50_CI = c(round(hdi(Baseline1_High,0.5)[[2]], digits = 2),
                              round(hdi(Baseline2_High,0.5)[[2]], digits = 2),
                        round(hdi(Baseline1_Low,0.5)[[2]], digits = 2),
                        round(hdi(Baseline2_Low,0.5)[[2]], digits = 2),
                        round(hdi(Acquisition1_High,0.5)[[2]], digits = 2),
                        round(hdi(Acquisition2_High,0.5)[[2]], digits = 2),
                        round(hdi(Acquisition1_Low,0.5)[[2]], digits = 2),
                        round(hdi(Acquisition2_Low,0.5)[[2]], digits = 2),
                        round(hdi(Extinction1_High,0.5)[[2]], digits = 2),
                        round(hdi(Extinction2_High,0.5)[[2]], digits = 2),
                        round(hdi(Extinction1_Low,0.5)[[2]], digits = 2),
                        round(hdi(Extinction2_Low,0.5)[[2]], digits = 2))


posterior_means =  posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")

posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
                                      "High Reward" = "High reward",
                                      "Low Reward" = "Low reward")

posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
                                "Acquisition1" = "Training1",
                                "Acquisition2" = "Training2",
                                "Baseline1" = "Baseline1",
                                "Baseline1" = "Baseline2",
                                "Extinction1" = "Test1",
                                "Extinction2" = "Test2")

# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$`dprime` = c(1,1,1,1,1,1)


# Raw data
# Prepare the dataset
data.plot = data.final

# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] = "Reward Phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] = "Reward Probability"

# rename conditions
data.plot$`Reward Phase` = recode(data.plot$`Reward Phase`,
                                "Acq1" = "Training1",
                                "Acq2" = "Training2",
                                "Bsln1" = "Baseline1",
                                "Bsln2" = "Baseline2",
                                "Ext1" = "Test1",
                                "Ext2" = "Test2")

data.plot$`Reward Probability` = recode(data.plot$`Reward Probability`,
                                      "High_Rew" = "High reward",
                                      "Low_Rew" = "Low reward")


# Violin + the model
# tiff("dprime_split.tiff", units="in", width=16, height=5, res=800)

f.sup.a = ggplot(data=data.plot, aes(x=`Reward Phase`, y=`dprime`, label=`Reward Probability`)) + 
  # violin of the raw data
  geom_violin(data=data.plot, color ="#636363",trim = TRUE,alpha="0.2") +
  # individual participants for the raw data
  # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
  geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
  # mean of the model
  stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`dprime`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
  # 95 credible intervals of the model
  geom_linerange(data = posterior_means,width=.1, aes(ymin=`dprime`-`dprime`+low_95_CI, ymax=`dprime`-`dprime`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
  # 50 credible intervals of the model
  geom_linerange(data = posterior_means,width=.1, aes(ymin=`dprime`-`dprime`+low_50_CI, ymax=`dprime`-`dprime`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
  scale_x_discrete(limits=c("Baseline1","Baseline2", "Training1","Training2", "Test1","Test2")) +
  facet_grid(.~`Reward Probability`) + 
  scale_y_continuous(limits = c(-1,3),breaks=seq(-1,3,0.5)) +
  labs(title="Sensitivity d′", y = "Sensitivity d′") + 
  theme(
    axis.line = element_line(size=1, colour = "black"),
    panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
    panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
    panel.border = element_blank(), 
    panel.background = element_blank(),
    plot.title = element_text(size = 18,  face = "bold",hjust = 0.5),
    text=element_text(colour="black", size = 14),
    axis.text.x=element_text(colour="black", size = 14),
    axis.text.y=element_text(colour="black", size = 14),
    axis.title=element_text(size=16,colour = "black")) + 
  theme(strip.background =element_rect(fill="white")) + 
  theme(strip.text = element_text(size = 16))

f.sup.a  

# dev.off()

Table of means across conditions

# Make a table with conditions
posterior_means = as.data.frame(c("Baseline1 High Reward", 
                                  "Baseline1 Low Reward", 
                                  "Baseline2 High Reward", 
                                  "Baseline2 Low Reward",
                                  "Acquisition1 High Reward", 
                                  "Acquisition1 Low Reward",
                                  "Acquisition2 High Reward", 
                                  "Acquisition2 Low Reward", 
                                  "Extinction1 High Reward", 
                                  "Extinction1 Low Reward",
                                  "Extinction2 High Reward",
                                  "Extinction2 Low Reward"))

names(posterior_means)[1] = "Condition"

posterior_means$Mean = c(paste(round(mean(Baseline1_High), digits = 2), " [", round(hdi(Baseline1_High)[[1]], digits = 2), " ", round(hdi(Baseline1_High)[[2]], digits = 2), "]"),
                         paste(round(mean(Baseline1_Low), digits = 2), " [", round(hdi(Baseline1_Low)[[1]], digits = 2), " ", round(hdi(Baseline1_Low)[[2]], digits = 2), "]"),
                         paste(round(mean(Baseline2_High), digits = 2), " [", round(hdi(Baseline2_High)[[1]], digits = 2), " ", round(hdi(Baseline2_High)[[2]], digits = 2), "]"),
                         paste(round(mean(Baseline2_Low), digits = 2), " [", round(hdi(Baseline2_Low)[[1]], digits = 2), " ", round(hdi(Baseline2_Low)[[2]], digits = 2), "]"),
                         paste(round(mean(Acquisition1_High), digits = 2), " [", round(hdi(Acquisition1_High)[[1]], digits = 2), " ", round(hdi(Acquisition1_High)[[2]], digits = 2), "]"),
                         paste(round(mean(Acquisition1_Low), digits = 2), " [", round(hdi(Acquisition1_Low)[[1]], digits = 2), " ", round(hdi(Acquisition1_Low)[[2]], digits = 2), "]"),
                         paste(round(mean(Acquisition2_High), digits = 2), " [", round(hdi(Acquisition2_High)[[1]], digits = 2), " ", round(hdi(Acquisition2_High)[[2]], digits = 2), "]"),
                         paste(round(mean(Acquisition2_Low), digits = 2), " [", round(hdi(Acquisition2_Low)[[1]], digits = 2), " ", round(hdi(Acquisition2_Low)[[2]], digits = 2), "]"),
                         paste(round(mean(Extinction1_High), digits = 2), " [", round(hdi(Extinction1_High)[[1]], digits = 2), " ", round(hdi(Extinction1_High)[[2]], digits = 2), "]"),
                         paste(round(mean(Extinction1_Low), digits = 2), " [", round(hdi(Extinction1_Low)[[1]], digits = 2), " ", round(hdi(Extinction1_Low)[[2]], digits = 2), "]"),
                         paste(round(mean(Extinction2_High), digits = 2), " [", round(hdi(Extinction2_High)[[1]], digits = 2), " ", round(hdi(Extinction2_High)[[2]], digits = 2), "]"),
                         paste(round(mean(Extinction2_Low), digits = 2), " [", round(hdi(Extinction2_Low)[[1]], digits = 2), " ", round(hdi(Extinction2_Low)[[2]], digits = 2), "]"))

names(posterior_means)[2] = "Mean [HDI]"

posterior_means =  posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")
                        
kable(posterior_means, caption = "Means per condition")
Means per condition
Reward Phase Reward Probability Mean [HDI]
Baseline1 High Reward 1.48 [ 1.24 1.71 ]
Baseline1 Low Reward 1.32 [ 1.09 1.53 ]
Baseline2 High Reward 1.6 [ 1.38 1.84 ]
Baseline2 Low Reward 1.47 [ 1.25 1.69 ]
Acquisition1 High Reward 1.54 [ 1.3 1.78 ]
Acquisition1 Low Reward 1.59 [ 1.37 1.81 ]
Acquisition2 High Reward 1.6 [ 1.35 1.83 ]
Acquisition2 Low Reward 1.48 [ 1.26 1.7 ]
Extinction1 High Reward 1.49 [ 1.24 1.72 ]
Extinction1 Low Reward 1.5 [ 1.28 1.71 ]
Extinction2 High Reward 1.49 [ 1.25 1.74 ]
Extinction2 Low Reward 1.55 [ 1.33 1.76 ]

Inference about the best model

Check the difference between baseline1 vs. baseline2 in high reward

Diff_Bsln1_Bsln2_High = Baseline2_High - Baseline1_High
plotPost(Diff_Bsln1_Bsln2_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

Check the difference between baseline1 vs. baseline2 in low reward

Diff_Bsln1_Bsln2_Low = Baseline2_Low - Baseline1_Low
plotPost(Diff_Bsln1_Bsln2_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

Check the difference between baseline2 vs. acquisition1 in high reward

Diff_Bsln2_Acq1_High = Baseline2_High - Acquisition1_High
plotPost(Diff_Bsln2_Acq1_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

Check the difference between baseline2 vs. acquisition1 in low reward

Diff_Bsln2_Acq1_Low = Acquisition1_Low - Baseline2_Low
plotPost(Diff_Bsln2_Acq1_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

Reaction times

Means - raw data

summary = ddply(data.final,.(ExpPhase,Condition),plyr::summarize,Mean=c(paste(round(mean(Hits.RTs, na.rm = TRUE), digits = 2), " [", round(hdi(Hits.RTs)[[1]], digits = 2), " ", round(hdi(Hits.RTs)[[2]], digits = 2), "]")))

names(summary) = c("Reward phase", "Reward probability", "Reaction times")

# rename conditions
summary$`Reward phase` = recode(summary$`Reward phase`,
                                "Acq1" = "Acquisition1",
                                "Acq2" = "Acquisition2",
                                "Bsln1" = "Baseline1",
                                "Bsln2" = "Baseline2",
                                "Ext1" = "Extinction1",
                                "Ext2" = "Extinction2")

summary$`Reward probability` = recode(summary$`Reward probability`,
                           "High_Rew" = "High",
                           "Low_Rew" = "Low")

summary = as.data.frame(summary)

kable(summary, caption = "Reaction times per condition")
Reaction times per condition
Reward phase Reward probability Reaction times
Baseline1 High 548.84 [ 479.43 613.76 ]
Baseline1 Low 548.43 [ 458.26 610.63 ]
Baseline2 High 544.34 [ 454.56 620.36 ]
Baseline2 Low 554.01 [ 479.48 632.8 ]
Acquisition1 High 521.4 [ 437.9 587.57 ]
Acquisition1 Low 542.34 [ 463.65 593.47 ]
Acquisition2 High 528.74 [ 462 598.58 ]
Acquisition2 Low 533.94 [ 479.38 618.25 ]
Extinction1 High 528.58 [ 457.88 596.17 ]
Extinction1 Low 536.54 [ 444.89 621 ]
Extinction2 High 529.3 [ 448.24 606 ]
Extinction2 Low 543.01 [ 450.11 617.44 ]

Plot

# Prepare the dataset
data.plot = data.final

# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] = "Reward phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] = "Reward probability"

# rename conditions
data.plot$`Reward phase` = recode(data.plot$`Reward phase`,
                                "Acq" = "Acquisition",
                                "Bsln" = "Baseline",
                                "Ext" = "Extinction")

data.plot$`Reward probability` = recode(data.plot$`Reward probability`,
                                      "High_Rew" = "High",
                                      "Low_Rew" = "Low")

# Pirate plot
  pirateplot(formula = Hits.RTs ~ `Reward phase` + `Reward probability`, # dependent~independent variables
             data=data.plot, # data frame
             main="Reaction times", # main title
             ylim=c(400,700), # y-axis: limits
             ylab="Reaction time", # y-axis: label
             theme=0, # preset theme (0: use your own)
             point.col="black", # points: color
             point.o=.3, # points: opacity (0-1)
             avg.line.fun = median,
             avg.line.col="black", # average line: color
             avg.line.lwd=2, # average line: line width
             avg.line.o=1, # average line: opacity (0-1)
             bean.b.col="black", # bean border, color
             bean.lwd=0.6, # bean border, line width
             bean.lty=1, # bean border, line type (1: solid; 2:dashed; 3: dotted; ...)
             bean.b.o=0.3, # bean border, opacity (0-1)
             bean.f.col="gray", # bean filling, color
             bean.f.o=.1, # bean filling, opacity (0-1)
             cap.beans=FALSE, # max and min values of bean densities are capped at the limits found in the data
             gl.col="gray", # gridlines: color
             gl.lty=2, # gridlines: line type (1: solid; 2:dashed; 3: dotted; ...)
             cex.lab=1, # axis labels: size
             cex.axis=1, # axis numbers: size
             cex.names = 1,
             bty="l", # plot box type
             back.col="white") # background, color

Statistics

# Set the working directory in order to load the models
# Set working directory
setwd(here("brms_models"))

model.full.RT.split = readRDS("model.full.RT.split.rds")
compare.waic.RT.split = readRDS("compare.RT.split.waic")
compare.RT.split.waic.weights = readRDS("compare.RT.split.waic.weights")
bR2.null.RT.split = readRDS("bR2.null.RT.split")
bR2.expphase.RT.split = readRDS("bR2.expphase.RT.split")
bR2.full.RT.split = readRDS("bR2.full.RT.split")

Model comparison with WAIC

print(compare.waic.RT.split)
## Output of model 'model.null.RT':
## 
## Computed from 12000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic  -2538.8 24.5
## p_waic        40.4  4.6
## waic        5077.7 49.0
## 
## Output of model 'model.expphase.RT':
## 
## Computed from 12000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic  -2527.6 26.3
## p_waic        80.2 10.5
## waic        5055.3 52.5
## 
## Output of model 'model.full.RT':
## 
## Computed from 12000 by 516 log-likelihood matrix
## 
##           Estimate   SE
## elpd_waic  -2426.8 24.5
## p_waic       137.6 13.2
## waic        4853.6 49.0
## 
## Model comparisons:
##                   elpd_diff se_diff
## model.full.RT        0.0       0.0 
## model.expphase.RT -100.8      13.0 
## model.null.RT     -112.0      13.5

Bayesian R squared

Null model

print(bR2.null.RT.split)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.4599668 0.02556954 0.4063619 0.5072804

Experiment phase model

print(bR2.expphase.RT.split)
##     Estimate Est.Error      Q2.5     Q97.5
## R2 0.5357306 0.0277915 0.4777634 0.5850535

Full model

print(bR2.full.RT.split)
##     Estimate  Est.Error      Q2.5     Q97.5
## R2 0.7298469 0.01868379 0.6897955 0.7633356

Checking the best model

Plotting the chains

# Plot chains
plot(model.full.RT.split, pars = "^b_", ask = FALSE, N=6)

Summary of the best model

# Summary of the best model
print(tidy(model.full.RT.split, par_type = "non-varying"), digits = 1)
##                              term estimate std.error lower upper
## 1                       Intercept    548.9         6   538   560
## 2                   ExpPhaseBsln2     -4.5         5   -13     4
## 3                    ExpPhaseAcq1    -27.4         5   -36   -19
## 4                    ExpPhaseAcq2    -20.1         5   -29   -11
## 5                    ExpPhaseExt1    -20.2         6   -30   -11
## 6                    ExpPhaseExt2    -19.5         6   -29   -10
## 7                ConditionLow_Rew     -0.5         7   -12    11
## 8  ExpPhaseBsln2:ConditionLow_Rew     10.1         7    -1    21
## 9   ExpPhaseAcq1:ConditionLow_Rew     21.4         7    10    33
## 10  ExpPhaseAcq2:ConditionLow_Rew      5.7         7    -6    17
## 11  ExpPhaseExt1:ConditionLow_Rew      8.4         7    -3    20
## 12  ExpPhaseExt2:ConditionLow_Rew     14.2         7     3    25

Plotting the best model

Sample from the posterior

# Analyzing the posterior and differences between conditions

post = posterior_samples(model.full.RT.split, "^b")


################################################ Baseline 1 ####

######### High reward
Baseline1_High = post[["b_Intercept"]]
######### Low reward
Baseline1_Low = post[["b_Intercept"]] + 
  post[["b_ConditionLow_Rew"]] 

################################################ Baseline 2 ####

######### High reward
Baseline2_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseBsln2"]] 
######### Low reward
Baseline2_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseBsln2"]] +
  post[["b_ConditionLow_Rew"]] +
  post[["b_ExpPhaseBsln2:ConditionLow_Rew"]]

################################################ Acquistion 1

######### High reward
Acquisition1_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq1"]] 
######### Low reward
Acquisition1_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq1"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ExpPhaseAcq1:ConditionLow_Rew"]]

################################################ Acquistion 2

######### High reward
Acquisition2_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq2"]] 
######### Low reward
Acquisition2_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseAcq2"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ExpPhaseAcq2:ConditionLow_Rew"]]

################################################ Extinction 1

######### High reward
Extinction1_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt1"]] 
######### Low reward
Extinction1_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt1"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ExpPhaseExt1:ConditionLow_Rew"]]

################################################ Extinction 2

######### High reward
Extinction2_High = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt2"]] 
######### Low reward
Extinction2_Low = post[["b_Intercept"]] + 
  post[["b_ExpPhaseExt2"]] + 
  post[["b_ConditionLow_Rew"]] + 
  post[["b_ExpPhaseExt2:ConditionLow_Rew"]]
# Data from the model

# make a data frame
posterior_conditions = melt(data.frame(Baseline1_High, Baseline2_High, Baseline1_Low, Baseline2_Low, Acquisition1_High, Acquisition1_Low, Acquisition2_High, Acquisition2_Low, Extinction1_High, Extinction1_Low, Extinction2_High, Extinction2_Low))

posterior_conditions =  posterior_conditions %>% separate(variable, c("Reward Phase", "Reward Probability"), "_", extra = "merge")

names(posterior_conditions)[3] = "Reaction time"

posterior_conditions$`Reward Probability` = recode(posterior_conditions$`Reward Probability`,
                                      "High" = "High reward",
                                      "Low" = "Low reward")

posterior_conditions$`Reward Phase` = recode(posterior_conditions$`Reward Phase`,
                                "Acquisition1" = "Training1",
                                "Acquisition2" = "Training2",
                                "Baseline1" = "Baseline1",
                                "Baseline1" = "Baseline2",
                                "Extinction1" = "Test1",
                                "Extinction2" = "Test2")

# Make a table with credible intervals
posterior_means = as.data.frame(c("Baseline1 High Reward", 
                                  "Baseline2 High Reward", 
                                  "Baseline1 Low Reward", 
                                  "Baseline2 Low Reward",
                                  "Acquisition1 High Reward", 
                                  "Acquisition2 High Reward",
                                  "Acquisition1 Low Reward", 
                                  "Acquisition2 Low Reward", 
                                  "Extinction1 High Reward", 
                                  "Extinction2 High Reward",
                                  "Extinction1 Low Reward",
                                  "Extinction2 Low Reward"))
names(posterior_means)[1] = "Condition"

posterior_means$low_95_CI = c(round(hdi(Baseline1_High)[[1]], digits = 2),
                              round(hdi(Baseline2_High)[[1]], digits = 2),
                        round(hdi(Baseline1_Low)[[1]], digits = 2),
                        round(hdi(Baseline2_Low)[[1]], digits = 2),
                        round(hdi(Acquisition1_High)[[1]], digits = 2),
                        round(hdi(Acquisition2_High)[[1]], digits = 2),
                        round(hdi(Acquisition1_Low)[[1]], digits = 2),
                        round(hdi(Acquisition2_Low)[[1]], digits = 2),
                        round(hdi(Extinction1_High)[[1]], digits = 2),
                        round(hdi(Extinction2_High)[[1]], digits = 2),
                        round(hdi(Extinction1_Low)[[1]], digits = 2),
                        round(hdi(Extinction2_Low)[[1]], digits = 2))

posterior_means$high_95_CI = c(round(hdi(Baseline1_High)[[2]], digits = 2),
                              round(hdi(Baseline2_High)[[2]], digits = 2),
                        round(hdi(Baseline1_Low)[[2]], digits = 2),
                        round(hdi(Baseline2_Low)[[2]], digits = 2),
                        round(hdi(Acquisition1_High)[[2]], digits = 2),
                        round(hdi(Acquisition2_High)[[2]], digits = 2),
                        round(hdi(Acquisition1_Low)[[2]], digits = 2),
                        round(hdi(Acquisition2_Low)[[2]], digits = 2),
                        round(hdi(Extinction1_High)[[2]], digits = 2),
                        round(hdi(Extinction2_High)[[2]], digits = 2),
                        round(hdi(Extinction1_Low)[[2]], digits = 2),
                        round(hdi(Extinction2_Low)[[2]], digits = 2))

posterior_means$low_50_CI = c(round(hdi(Baseline1_High,0.5)[[1]], digits = 2),
                              round(hdi(Baseline2_High,0.5)[[1]], digits = 2),
                        round(hdi(Baseline1_Low,0.5)[[1]], digits = 2),
                        round(hdi(Baseline2_Low,0.5)[[1]], digits = 2),
                        round(hdi(Acquisition1_High,0.5)[[1]], digits = 2),
                        round(hdi(Acquisition2_High,0.5)[[1]], digits = 2),
                        round(hdi(Acquisition1_Low,0.5)[[1]], digits = 2),
                        round(hdi(Acquisition2_Low,0.5)[[1]], digits = 2),
                        round(hdi(Extinction1_High,0.5)[[1]], digits = 2),
                        round(hdi(Extinction2_High,0.5)[[1]], digits = 2),
                        round(hdi(Extinction1_Low,0.5)[[1]], digits = 2),
                        round(hdi(Extinction2_Low,0.5)[[1]], digits = 2))

posterior_means$high_50_CI = c(round(hdi(Baseline1_High,0.5)[[2]], digits = 2),
                              round(hdi(Baseline2_High,0.5)[[2]], digits = 2),
                        round(hdi(Baseline1_Low,0.5)[[2]], digits = 2),
                        round(hdi(Baseline2_Low,0.5)[[2]], digits = 2),
                        round(hdi(Acquisition1_High,0.5)[[2]], digits = 2),
                        round(hdi(Acquisition2_High,0.5)[[2]], digits = 2),
                        round(hdi(Acquisition1_Low,0.5)[[2]], digits = 2),
                        round(hdi(Acquisition2_Low,0.5)[[2]], digits = 2),
                        round(hdi(Extinction1_High,0.5)[[2]], digits = 2),
                        round(hdi(Extinction2_High,0.5)[[2]], digits = 2),
                        round(hdi(Extinction1_Low,0.5)[[2]], digits = 2),
                        round(hdi(Extinction2_Low,0.5)[[2]], digits = 2))


posterior_means =  posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")

posterior_means$`Reward Probability` = recode(posterior_means$`Reward Probability`,
                                      "High Reward" = "High reward",
                                      "Low Reward" = "Low reward")

posterior_means$`Reward Phase` = recode(posterior_means$`Reward Phase`,
                                "Acquisition1" = "Training1",
                                "Acquisition2" = "Training2",
                                "Baseline1" = "Baseline1",
                                "Baseline1" = "Baseline2",
                                "Extinction1" = "Test1",
                                "Extinction2" = "Test2")

# Create this variable in order to have it for plotting (this is a quick hack: in plotting the CIs this variable is added and then subtracted)
posterior_means$`Reaction time` = c(555,555,555,555,555,555)


# Raw data
# Prepare the dataset
data.plot = data.final

# rename variables
colnames(data.plot)[colnames(data.plot)=="ExpPhase"] = "Reward Phase"
colnames(data.plot)[colnames(data.plot)=="Condition"] = "Reward Probability"

# rename conditions
data.plot$`Reward Phase` = recode(data.plot$`Reward Phase`,
                                "Acq1" = "Training1",
                                "Acq2" = "Training2",
                                "Bsln1" = "Baseline1",
                                "Bsln2" = "Baseline2",
                                "Ext1" = "Test1",
                                "Ext2" = "Test2")

data.plot$`Reward Probability` = recode(data.plot$`Reward Probability`,
                                      "High_Rew" = "High reward",
                                      "Low_Rew" = "Low reward")

data.plot$`Reaction time` = data.plot$Hits.RTs


# Violin + the model
# tiff("RTs_split.tiff", units="in", width=16, height=5, res=800)

f.sup.b = ggplot(data=data.plot, aes(x=`Reward Phase`, y=`Reaction time`, label=`Reward Probability`)) + 
  # violin of the raw data
  geom_violin(data=data.plot, color ="#636363",trim = TRUE,alpha="0.2") +
  # individual participants for the raw data
  # geom_point(binaxis='y', stackdir='center', dotsize=0.5, alpha="0.2") + 
  geom_jitter(dotsize=0.9, alpha="0.6",width = 0.1,color="#bdbdbd")+
  # mean of the model
  stat_summary(data = posterior_conditions, aes(x=`Reward Phase`, y=`Reaction time`, group=`Reward Probability`),fun.y = mean, geom="line",size=0.8, colour = "#636363")+ 
  # 95 credible intervals of the model
  geom_linerange(data = posterior_means,width=.1, aes(ymin=`Reaction time`-`Reaction time`+low_95_CI, ymax=`Reaction time`-`Reaction time`+high_95_CI), size = 4, color = "#74a9cf", alpha = 1) +
  # 50 credible intervals of the model
  geom_linerange(data = posterior_means,width=.1, aes(ymin=`Reaction time`-`Reaction time`+low_50_CI, ymax=`Reaction time`-`Reaction time`+high_50_CI), size = 4, color = "#045a8d", alpha = 1) +
  scale_x_discrete(limits=c("Baseline1","Baseline2", "Training1","Training2", "Test1","Test2")) +
  facet_grid(.~`Reward Probability`) + 
  scale_y_continuous(limits =c(400,700),breaks=seq(400,700,50)) +
  labs(title="Reaction time", y = "Reaction time (ms)") + 
  theme(
    axis.line = element_line(size=1, colour = "black"),
    panel.grid.major = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
    panel.grid.minor = element_line(colour = "grey",size = 0.1,linetype = "dashed"),
    panel.border = element_blank(), 
    panel.background = element_blank(),
    plot.title = element_text(size = 18,  face = "bold",hjust = 0.5),
    text=element_text(colour="black", size = 14),
    axis.text.x=element_text(colour="black", size = 14),
    axis.text.y=element_text(colour="black", size = 14),
    axis.title=element_text(size=16,colour = "black")) + 
  theme(strip.background =element_rect(fill="white")) + 
  theme(strip.text = element_text(size = 16))
  
# dev.off()

f.sup.b

Table of means across conditions

# Make a table with conditions
posterior_means = as.data.frame(c("Baseline1 High Reward", 
                                  "Baseline1 Low Reward", 
                                  "Baseline2 High Reward", 
                                  "Baseline2 Low Reward",
                                  "Acquisition1 High Reward", 
                                  "Acquisition1 Low Reward",
                                  "Acquisition2 High Reward", 
                                  "Acquisition2 Low Reward", 
                                  "Extinction1 High Reward", 
                                  "Extinction1 Low Reward",
                                  "Extinction2 High Reward",
                                  "Extinction2 Low Reward"))

names(posterior_means)[1] = "Condition"

posterior_means$Mean = c(paste(round(mean(Baseline1_High), digits = 2), " [", round(hdi(Baseline1_High)[[1]], digits = 2), " ", round(hdi(Baseline1_High)[[2]], digits = 2), "]"),
                         paste(round(mean(Baseline1_Low), digits = 2), " [", round(hdi(Baseline1_Low)[[1]], digits = 2), " ", round(hdi(Baseline1_Low)[[2]], digits = 2), "]"),
                         paste(round(mean(Baseline2_High), digits = 2), " [", round(hdi(Baseline2_High)[[1]], digits = 2), " ", round(hdi(Baseline2_High)[[2]], digits = 2), "]"),
                         paste(round(mean(Baseline2_Low), digits = 2), " [", round(hdi(Baseline2_Low)[[1]], digits = 2), " ", round(hdi(Baseline2_Low)[[2]], digits = 2), "]"),
                         paste(round(mean(Acquisition1_High), digits = 2), " [", round(hdi(Acquisition1_High)[[1]], digits = 2), " ", round(hdi(Acquisition1_High)[[2]], digits = 2), "]"),
                         paste(round(mean(Acquisition1_Low), digits = 2), " [", round(hdi(Acquisition1_Low)[[1]], digits = 2), " ", round(hdi(Acquisition1_Low)[[2]], digits = 2), "]"),
                         paste(round(mean(Acquisition2_High), digits = 2), " [", round(hdi(Acquisition2_High)[[1]], digits = 2), " ", round(hdi(Acquisition2_High)[[2]], digits = 2), "]"),
                         paste(round(mean(Acquisition2_Low), digits = 2), " [", round(hdi(Acquisition2_Low)[[1]], digits = 2), " ", round(hdi(Acquisition2_Low)[[2]], digits = 2), "]"),
                         paste(round(mean(Extinction1_High), digits = 2), " [", round(hdi(Extinction1_High)[[1]], digits = 2), " ", round(hdi(Extinction1_High)[[2]], digits = 2), "]"),
                         paste(round(mean(Extinction1_Low), digits = 2), " [", round(hdi(Extinction1_Low)[[1]], digits = 2), " ", round(hdi(Extinction1_Low)[[2]], digits = 2), "]"),
                         paste(round(mean(Extinction2_High), digits = 2), " [", round(hdi(Extinction2_High)[[1]], digits = 2), " ", round(hdi(Extinction2_High)[[2]], digits = 2), "]"),
                         paste(round(mean(Extinction2_Low), digits = 2), " [", round(hdi(Extinction2_Low)[[1]], digits = 2), " ", round(hdi(Extinction2_Low)[[2]], digits = 2), "]"))

names(posterior_means)[2] = "Mean [HDI]"

posterior_means =  posterior_means %>% separate(Condition, c("Reward Phase", "Reward Probability"), " ", extra = "merge")

kable(posterior_means, caption = "Means per condition")
Means per condition
Reward Phase Reward Probability Mean [HDI]
Baseline1 High Reward 548.86 [ 535.97 561.35 ]
Baseline1 Low Reward 548.38 [ 535.83 560.97 ]
Baseline2 High Reward 544.34 [ 531.22 558.49 ]
Baseline2 Low Reward 553.98 [ 540.67 567.69 ]
Acquisition1 High Reward 521.42 [ 508.48 533.66 ]
Acquisition1 Low Reward 542.35 [ 530.05 555.45 ]
Acquisition2 High Reward 528.74 [ 515.92 541.36 ]
Acquisition2 Low Reward 533.91 [ 521.41 547.24 ]
Extinction1 High Reward 528.64 [ 514.39 542.24 ]
Extinction1 Low Reward 536.51 [ 520.49 551.37 ]
Extinction2 High Reward 529.32 [ 516.53 543.7 ]
Extinction2 Low Reward 543.01 [ 528.56 557.28 ]

Inference about the best model

Check the difference between baseline1 vs. baseline2 in high reward

Diff_Bsln1_Bsln2_High = Baseline2_High - Baseline1_High
plotPost(Diff_Bsln1_Bsln2_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

Check the difference between baseline1 vs. baseline2 in low reward

Diff_Bsln1_Bsln2_Low = Baseline2_Low - Baseline1_Low
plotPost(Diff_Bsln1_Bsln2_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

Check the difference between baseline2 vs. acquisition1 in high reward

Diff_Bsln2_Acq1_High = Baseline2_High - Acquisition1_High
plotPost(Diff_Bsln2_Acq1_High, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

Check the difference between baseline2 vs. acquisition1 in low reward

Diff_Bsln2_Acq1_Low = Acquisition1_Low - Baseline2_Low
plotPost(Diff_Bsln2_Acq1_Low, xlab = "", col = "#b3cde0", showCurve = FALSE, cex = 1, compVal = 0)

Creating figures for the manuscript

library(cowplot)
# Figure 2
tiff("C:/Users/igrahek/Documents/Studies/SSVEP Reward - Soren & Antonio/Experiment 1/SSVEP_reward/figures/Figure2.tiff", units="in", width=15, height=10, res=300)

ggdraw() +
    draw_plot(f.2.a, x = 0.03, y = 0.5, width = .45, height = .45) +
    draw_plot(f.2.b, x = 0.53, y = .5, width = .45, height = .45) +
    draw_plot(f.2.c, x = 0.03, y = 0, width = 0.45, height = 0.45) +
    draw_plot(f.2.d, x = 0.53, y = 0, width = 0.45, height = 0.45) +

  draw_plot_label(label = c("A", "B", "C","D"), size = 20,
                  x = c(0, 0.5, 0,0.5), y = c(1, 1, 0.5,0.5))

dev.off()
## png 
##   2
# Figure 1 Supplementary materials
tiff("C:/Users/igrahek/Documents/Studies/SSVEP Reward - Soren & Antonio/Experiment 1/SSVEP_reward/figures/FigureSup1.tiff", units="in", width=15, height=10, res=300)

ggdraw() +
    draw_plot(f.sup.a, x = 0.03, y = 0.5, width = .9, height = .45) +
    draw_plot(f.sup.b, x = 0.03, y = 0, width = .9, height = .45) +

  draw_plot_label(label = c("A", "B"), size = 20,
                  x = c(0, 0), y = c(1,0.5))

dev.off()
## png 
##   2